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Nomenclature 


c. 

1 


s coefficients in polynomial or nondimensional prediction functions 


* .th 

D = value of dependent variable (impact damage) at i data point 


J MAX 


maximum diameter of the bumper hole 

D = minimum diameter of the bumper hole 
MIN 

n =s average diameter of the hole in the MLI 
MLI 

D * diameter of projectile 
P 

D = average diameter of the pressure wall hole, 
pw 

D = bumper stand-off distance 
s 

D = estimated value of dependent variable 

interpolation or prediction point 

E = elastic modulus of the bumper plate material 

M = number of data points in database or number of material properties in 

each record of the materials data file 

N = number of independent variables (impact parameters) 

rn# = random number used to generate test function 

„2 


(impact damage) at an 


R. 

i 


pw 


\ • « 

J.i 


= coefficient of determination 

= distance from data point to interpolation or prediction point 
= length of influence of a data point 
» bumper thickness 
= pressure wall thickness 

= speed of sound in the bumper material = J E /p 
= j** 1 coordinate (independent variable) of i** 1 data point 

j t ^ 1 coordinate (independent variable) of interpolation 


X j,INT = 
prediction point 

Ax = x - x., T ._ 
j’i j’i j INT 

V = projectile velocity 

0 = impact angle 

p = mass density of the bumper plate material 

0 = weighting factor of a data point 


or 
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1. Introduction 

There are many engineering applications where predictions of the 
behavior of a physical system must be made based on a database of 

experimental results. In these instances, either the phenomenon is too 

complicated to treat analytically or numerically, or the funding, expertise 
or time required to do so is not available. Empirical approaches of this 
nature have always played a fundamental role in engineering design. 

The purpose of this project was to develop a DOS microcomputer-based 
computer program to empirically predict hypervelocity impact damage to the 
Space Station from space debris. The main goal was to predict damage to the 
multilayer insulation (MLI). However, to extend the usefulness of the 

program, damage to other components of the Space Station wall are predicted 
as well. The program is intended to be an easy to use design tool for trade 

studies on debris protection strategies for the Space Station. The 

predictions are made based on a database of experimental results. 

MSFC has a light gas gun that can launch 2.5-12.7 mm projectiles at 
speeds of 2-8 km/s. 1 Work is currently in progress at MSFC to qualify the 
orbital debris protection system under development by Boeing for Space 
Station Freedom. A schematic of the protection system is shown in Fig. 1.1. 
It is based on the classical sacrificial bumper approach first suggested by 
Whipple. The purpose of the bumper is to break-up or ideally vaporize the 
projectile (space debris or micrometeoroids) so that the pressurized 

spacecraft behind the bumper is impacted with a series of fine particles 
rather than a single large particle. 

The parameters associated with the impact data are illustrated in Fig. 
1.1. The projectiles were initially spherical and typically constructed of 
1100 aluminum. The bumper and the pressure walls were typically made from 
6061-T6 and 2219-T87 aluminum, respectively. Some tests have been run 
with different materials. If a number of different materials were used for 
the bumper and the other components in the same database, then the number of 
independent degrees of freedom would be increased dramatically. This is 
because material properties, such as densities and melting points, would also 
have had to be accounted for. 

There are five computer programs that were developed for this project. 
Details of how to use the program are provided in section 2. A main program 
called MLIBLAST serves as a shell to run the other four programs. A program 
called DATABASE is provided to assist the user in creating and adding to the 
database file of experimental results. 
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Schematic drawing of impact specimen. 


The remaining three programs provide predictions of impact damage to the 
bumper, the multilayer insulation (MLI), and the pressure wall plate. Program 
INVRMETH uses a unique prediction technique, called the inverse R method, 
that was developed for the purposes of this project. The theoretical basis of 
this method is described in section 3. 

As described in section 4, program POLYMETH makes predictions by fitting 
simple polynomials through a subset of the data points. A more sophisticated 
form of polynomial prediction technique using the isoparametric f ormulation 
of the finite element method (FEM) as a basis was also attempted during the 
course of this project. This FEM based software was found to be somewhat 
unreliable for making predictions from the impact data that is currently 
available and so it was not included with this software. The interested 
reader can consult Appendix 1 for a discussion of this method. 

The last prediction program, NONDIMEN makes predictions based on 
nondinensionalized functions that were developed by others and extended by 
the author f or application here. These f unctions are described in section 5. 
The relative accuracies of these three prediction schemes are compared using 
an actual impact data set in section 6. 

Lists of conclusions and recommendations derived from this research 
project are given in section 7. 

A listing of the Microsoft BASIC source code is provided in Appendix 2 
or is provided on the computer disks. 
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2. Software User Guide 

The software developed for this project was written in Microsoft BASIC 
for DOS and compiled using the Microsoft BASIC Professional Development 
System compiler (version 7.0) for DOS based micro-computers. Approximately 
0.5 MB of hard disk space, am EGA or VGA graphics card and monitor, and an 
Intel 80286, 80386 or 80486 CPU is required to run the software. A math 
coprocessor is desirable, but not required. The source code may be modified 
amd recompiled using Microsoft QuickBASIC if desired. 

The software is provided on two 5.25", 360K computer disks. An annotated 
listing of the contents of the computer disks follows: 

DISK 1 

DATABASE. BAS - source code for the data base prograun (ASCII). 

DATABASE.EXE - compiled version of the database program. 

INVRMETH.BAS - source code for the inverse R method damage prediction program 
(ASCII). 

INVRMETH.EXE - compiled version of the inverse R method damage prediction 
program. 

MATERIAL.DAT - a typical database file of material properties which is used 
by the INVRMETH program (ASCII). 

MLI.DAT - a typical database file of experimental results (ASCII). 

MLIBLAST.BAS - source code for the main program that runs the other programs 
(ASCII). 

MLIBLAST.EXE - compiled version of the main program. 

DISK 2 

NONDIMEN. BAS - source code for the nondimensional functions damage prediction 
program (ASCII). 

NONDIMEN. EXE - compiled version of the nondimensional functions damage 
prediction program (ASCII). 

POLYMETH.BAS - source code for the polynomial functions damage prediction 
program (ASCII). 

POLYMETH.EXE - compiled version of the polynomial functions damage prediction 
program. 

The software is installed by first creating a sub directory on the hard 
disk and then copying all of the files from the two computer disks supplied 
into that subdirectory. If disk space is a problem then the source code files 
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(filename. BAS,) need not be copied. The program is started by typing MLIBLAST 
and following the prompts. More details on the program prompts are given 
below, but first the database files MATERIAL.DAT and MLI.DAT will be 


discussed. 

MATERIAL.DAT is an example of a typical materials data file. Any valid 
DOS name can be used for this file. Thus, the user may have several of this 
type of data file in a directory for different purposes. A file of this 
nature is required while running the inverse R program. The materials data 
file is an ASCII file that can be created and modified using any standard 
text editor. The format of the file is as follows: 


material property 1 name string 
material property 2 name string 


LISTING OF NAMES OF MATERIAL 
► PROPERTIES TO BE MODELED (MAXIMUM OF 10) 
(25 CHARACTERS MAX) 


material property M name string 

{ 

material 1 name string 

material property 1 for material 1 

material property 2 for material 1 

► TYPICAL DATA RECORD 


material property M for material 1 

} 

► ANY NUMBER OF DATA RECORDS MAY BE USED 


A typical material data file is provided on the computer disks and is called 

MATERIAL.DAT. This file is reproduced below: 

Density (lb/in~3) 

Elastic Mod. (lb/in*2) 

Ultimate Strgth (lb/ in" 2) 

Sp. Heat (BTU/db-deg R)) 

Melting Temp (deg R) 

{ 

1100 

9.780E-2 
1.000E 7 
1.600E4 
2.140E-1 
1.680E3 
} 

{ 

2219-T87 

1.030E-1 

1.050E7 

6.300E4 

2.050E-1 

1.680E3 
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} 

{ 

6061-T6 

9.800E-2 

9.900E6 

4.200E4 

2.100E-1 

1.680E3 

} 

MATERIAL.DAT is set up to model the material properties: density, elastic 
modulus, ultimate strength, specific heat, and melting temperature. Other 
physical properties can be used to a maximum of 10. The units do not have to 
be included in the material property name string. MATERIAL.DAT contains three 
records of material data for materials: 1100, 2219-T87, and 6061-T6. The 

material names are treated as string variables and thus can be any 
combination of numbers and letters. Any number of records of material data 
may be included. The order of the material properties must be the same in 
every record and be as the material property name strings are listed. For 
instance, referring to file MATERIAL.DAT, the specific heat of material 
2219-T87 is 2.050E-1. 

The purpose of the material properties database file is to provide an 
efficient, yet very flexible scheme for inputting material property data into 
the inverse R method computer program. The user can easily change the 
material properties to be modeled without disturbing the database file of 
experimental results. If the materials used for the projectile, bumper and 
pressure wall do not vary in the database, then the contents of the material 
properties database file will have no effect on the damage predicted by the 
inverse R method program. The polynomial function method program assumes that 
the material properties do not vary in the database. The nondimensional 
function method program assumes that the material properties of the 
projectile and pressure wall do not vary in the database and inputs material 
properties associated with the bumper directly. 

The other database file required for running the programs of MLIBLAST is 
associated with the experimental data. This file can be created (and 
enlarged) by running the DATABASE program from inside MLIBLAST or it can be 
created "by hand" using any standard text editor since it is an ASCII file. 
This file can be given any valid DOS file name. Currently, up to 100 data 
records can be placed in this file. The format for this file is as follows: 

{ 

Test Number String 

Test Agency String (source of data) 

Test Date String 
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Bumper Material String (same format as in material database file) 
Bumper Thickness 
Bumper Stand-Off 

Pressure Wall Material (same format as in material database file) 
Pressure Wall Thickness 

Projectile Material (same format as in material database file) 

Projectile Diameter 

Impact Angle 

Pro jectile Velocity 

Bumper Hole Maximum Diameter (Ma jor Axis) Dimension 
Bumper Hole Minimum Diameter (Minor Axis) Dimension 
MLI Mean Hole Diameter 
MLI Mass Loss 

Pressure Wall Hole Maximum Diameter (Major Axis) Dimension 
Pressure Wall Hole Minimum Diameter (Minor Axis) Dimension 

) 


► AS MANY AS 99 MORE DATA RECORDS 


MLI. DAT is an example of an experimental database file. This file is provided 
on the computer disks. It contains information on the specimens recently used 
for thermal testing in Sunspot Thermal Vacuum Chamber. To help understand the 
format information given above, the first record of MLI. DAT is presented 
below for comparison: 

{ 

1012 

MSFC 

05/08/90 

6061-T6 

.08 

4 

2219-T87 

.125 

1100 

.313 

0 

6.72 

.729 

.729 

2.2 

.938 

.6 

.15 

} 


An overview of the menu choices available to the user of MLIBLAST will 
now be discussed. The program is started by typing MLIBLAST. The user is then 
provided with three options: 

1. Add data to, or create a new experimental results database file. Selecting 
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this option will cause program DATABASE to rim. 


2. Make a prediction. This option involves running one of the three 
prediction programs: INVRMETH, NONDIMEN, or POLYMETH. 

3. Quit MLIBLAST 

The steps associated with running each of these programs will now be 
considered. 

PROGRAM DATABASE: 

Step 1 - Enter the name of an experimental results database file. Any valid 
DOS name can be used. If this file already exists, then the new data records 

will be appended to the end of it. MLI.DAT is an example of an experimental 

results data file. This file was provided on the computer disks. 

Step 2 - Enter the appropriate data at the prompts. Press ENTER after the 

data has been typed in. If you make a mistake, then press the F10 function 

key and then the ENTER key to redo the data input.- 

Step 3 - You will be prompted as to whether to write your previously entered 
data record information to your database file. This provides another way of 
not saving an input data record with errors. You will also be prompted as to 
whether to enter another data record. A response of n will cause you to exit 
from the database program. Note - the database file created is an ASCII file 
which can be edited with a standard text editor. Additional data records can 
be added to the experimental database file using the text editor (instead of 
program DATABASE) if so desired. 

PROGRAM INVRMETH (SEE SECTIOH 3 FOR MORE DETAILS ON PROGRAM INVRMETH) 

Step 1 - Input the names of the experimental database file and the material 
database file. The program will then read these files and present a summary 
of their contents on two computer screens. These summary screens are intended 
to help the user determine if the contents of the database file are 
appropriate for the desired prediction. 

Step 2 - Select the quantity for which a prediction is to be made. This 
program is designed to make predictions for: bumper hole maximum and minimum 
hole dimensions, MLI average hole diameter, MLI mass loss, and pressure wall 
maximum and minimum hole dimensions. 

Step 3 - Input the impact parameters (such as projectile diameter) associated 
with the desired prediction. Default values are provided in square brackets 
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for all inputs here except for impact angle. A default value is selected by 
simply pressing ENTER. The magnitude of the input impact parameter relative 
to the database average is indicated in round brackets. For instance if the 
projectile velocity for the prediction is twice that of the average 
projectile velocity in the experimental database file then the number 2 would 
appear in round brackets. Ideally, prediction parameter values should be 
close to the database average if reliable predictions are to be made. The 

round bracket numbers are intended to help the user assess the reliability of 
the prediction. 

Step 4 - Review the results of the prediction. Here, the value for the 

prediction is given and the location of the prediction point along the 
prediction vector (see section 3) is indicated. Information on the polynomial 
fit through the 10 interpolation points (see section 3) is also provided. The 
user can also review the results of the prediction graphically. Here, the 

variation of the function to be predicted along the prediction vector is 
illustrated to assist the user in assessing the reliability of the 

prediction. If the function being predicted varies in an erratic fashion 
along the prediction vector then the prediction may be unreliable. 

PROGRAM POLYMETH (SEE SECTION 4 for more details on program POLYMETH) 

Step 1 - A warning screen is displayed indicating that the program only 
models the parameters: bumper thickness, projectile diameter, projectile 

velocity, and impact angle. 

Step 2 - Input experimental results database filename. No material data file 
name is requested since it is assumed that material types will not vary in 
the database. The program will then read in the contents of the experimental 
results database file and display a summary of this data on the screen so 
that the user may assess its suitability with respect to the required 
predictions. 

Step 3 - Select the desired prediction such as MLI hole diameter. 

Step 4 - Input the impact parameters (such as bumper thickness) associated 
with the prediction. 

Step 5 - The program now attempts to fit a complete linear polynomial through 
random subsets of the data as described in section 4. If the experimental 
data does not "span" the impact parameter space very well then the program 
may take a long time or may not be able to to find a solution and make a 
prediction. In this case the user may press function key FI to quit or may 
try to fit incomplete polynomial functions. To attempt to fit a "simpler" 
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function (only terms to second order retained) the user should press function 
key F2. Pressing F3 will cause the program to fit the "simplest" f unction 
which retains only terms to first order. 

Step 6 - Analyze the results. The program will attempt to find five fits of 
the polynomial to random subsets of the experimental data. These will be 
displayed on the screen along with the average nondimensionalized distance 
(see section 4) of the data points used for the function fit. Intuitively, 
the prediction with the smallest average distance (data points closest to 
where the prediction is required) should be the most reliable. However, this 
may not be so if one of those data points happened to contain a large 
experimental error. Accordingly, a weighted average of the three predictions 
with the lowest mean distances is also presented. The weighting is based on 
the mean distances as described in section 4. 

PROGRAM NONDIMEN (SEE section 5 for more details on program nondimen) 

Step 1 - Input experimental database file name. No material database file 

name is required- because it is assumed that material types will not vary in 

the database. 

Step 2 - Input bumper elastic modulus (equal to 70E3 MPa for aluminum) and 

3 

input the mass density of the bumper (equal to 2710 kg/m for aluminum). 

After these values have been input, information summarizing the contents of 
the experimental database file will be shown on the screen. 

Step 3 - Parameters for the nonlinear function coefficient optimizer are 
input. The purpose of these parameters and recommended magnitudes are 
displayed on the computer screen. 

Step 4 - An this time an iterative procedure is invoked to adjust the 

prediction function coefficients such that the coefficient of determination 

2 2 
(R ) is maximized. During this process the R values are printed on the 

screen so that the user can assess the suitability of the functional form of 

prediction equations for fitting the experimental data. 

Step 5 - On completion of the optimization process, the function coefficients 
and R values are displayed on the screen to further assist the user in 
assessing the goodness of fit between the functions and the experimental 
data. 

Step 5 - Input prediction parameters (such as bumper thickness) and make 
predictions. 

In the next three sections more details on the prediction schemes are 
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presented. In section 6, the three prediction schemes are compared using the 
experimental data set associated with the impact specimens that were tested 
in the Sunspot Thermal Vacuum Chamber. 
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3. The Inverse R Prediction Technique 

The usual procedure for making predictions from experimental data is to 

assume some form for the equation relating the independent variables to the 

dependent variable. A function of this nature is described in section 5 of 

this report. The equation typically contains empirical coefficients, the 

3-8 

values of which are determined from a fit to the experimental data. The 
method of least squares (maximizing the coefficient of determination, R*2) is 
an example of a popular technique for obtaining the coefficients from the 
experimental data. The final result is a closed form equation for making 
predictions. 

This approach has been found to work very well for many engineering 
applications, however there are some disadvantages. A suitable form for the 
prediction equation must be developed. This is often difficult. Incorporating 
additional independent variables in an existing equation can pose problems. 
Usually, a well defined procedure for taking into account new experimental 
data is not put in place. Generally, a single set of empirical coefficients 
are used to make predictions over a fairly wide range of values of the 
independent variables. Thus, the best data in a database for making a 
prediction with a particular set of independent variables may not be used to 
best advantage. Also, it is usually difficult to assess the accuracy of the 

prediction. 

In this section, a new method (called inverse R method) for making 

empirical predictions based on experimental data is discussed. The method 

uses a very general form of prediction equation that can be applied in the 
same manner to all problems. Thus, the user is not required to develop a 
suitable form for the prediction equation and additional independent 
variables can be easily incorporated. The new method is designed to work off 
a database that can be continuously updated as new experimental data becomes 
available. The method automatically takes advantage of the most appropriate 
data in the database for a given set of independent variables. The method 

provides diagnostics for assessing the accuracy of the prediction. 

The new technique consists of four main steps which will now be 
described. 

Step 1. Normalize the Independent Variables. 

In general, the independent variables will vary greatly in magnitude. In 
hypervelocity impact work, dimensions can be of order 10 and velocities of 
order 10^. The new technique requires that all variables be of the same order 
of magnitude. This was accomplished by scaling the independent variables such 
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that their mean value is equal to unity. Other scaling methods could perhaps 

be used to improve the accuracy of this technique. For instance, the 

variables could be scaled such that predicted values of points in the 
database more closely match the measured values. This scaling technique was 
not tested. The dependent variables need not be scaled. 

This technique works off a database that can and should be kept updated 
with the latest experimental data. Thus, the scaling factors will change as 
time progresses and the size of the database increases. 

Step 2. Select a Series of Points in the Data Domain For Interpolation. 

Two general requirements for prediction schemes are: the method should 
be capable of smoothing the data to (hopefully) cancel out the random scatter 
typically present in experimental measurements and the technique should allow 
for making reliable predictions outside of the domain of the measured data. 
Here, these requirements are satisfied by using the data to make ten 

interpolations from within the domain of the data, which are then used for 

predicting the dependent variable at some point of interest. The ten 
"interpolation" points should provide for sufficient smoothing of the data 
and also capture the trend characteristics of the data for extrapolation 
purposes, if an extrapolation is required. The number of interpolation points 
to use was selected on the basis of trial and error. Note, in some cases 
extrapolation can produce misleading results regardless of the extrapolation 
technique used. 

Fig. 3.1 provides an illustration of how the interpolation points are 

selected for a hypothetical case with two independent variables. An identical 
approach is used for the case of an arbitrary number of independent 

variables. In Fig. 3.1, the independent variables are in the plane of the 
page and the dependent variable takes the form of a surface out of the plane 
of the page. 

First, a "prediction" vector is drawn from the origin through the point 

in the domain where a prediction of the dependent variable is required, which 
is called the "target" point. Then the "min" and "max" points (Fig. 3.1) are 
located on the prediction vector by considering the intersection points of 
perpendiculars from the data points to the prediction vector. The closest 
intersection point to the origin defines the min point, and that of the 
farthest, the max point. Ten equally spaced points (interpolation points) on 

the prediction vector between the min and max point are then used for the 
next step in the prediction process. If the target point lies between the min 
and max points then an interpolation is required, otherwise an extrapolation 
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is required. 

Step 3. Estimate Values of the Dependent Variable at Interpolation Points. 

Next, values for the dependent variable must be estimated at the ten 
interpolation points. This is done as indicated in the following equation: 


D = 



D. 

1 



1 

R^ 1 


l 


(3.1) 


The R. are determined by the usual formula for determining the "distance" 
between two points in a multidimensional space: 


N 

R^f = y (x. .-x. 

i L j,i j.INT 
j=l 


(3.2) 


The need for scaling the independent variables is evident from considering 
the form of Eq. (3.2). 

The form of Eq. (3.1) will now be considered. It is assumed that if all 
measured data points are the same "distance" R from an interpolation point 
then all the measured data should be given equal weight. This situation is 
illustrated for the case of two independent variables (N = 2) in Fig. 3.2. 
This can be interpreted as saying that each data point has some 
"characteristic length of influence", S, that subtends an angle 0 = S/R = 
S/R N_1 as indicated in Fig. 3.2. The 9 can be taken as the weighting factor. 
For the constant R case shown in Fig. 3.2, all data points would be given the 
same weight. Fig. 3.3 illustrates the case for which the data points are 
considered to be equally valid (same S), but are located different 

"distances" from the interpolation point. Here, the weighting factors will be 
of the form 0. = S/R^ _1 , and thus data points closer to the interpolation 
point will be given a higher weight. The value of the dependent variable at 
the interpolation point can be estimated from D = Z0.D./Z0. which leads to 
Eq. (3.1) and hence this technique is given the name Inverse R Method. Note 
that a value for S is not required as it cancels out of the equation. 

The three dimensional (three independent variables) application of this 
procedure leads to equations identical in form to those used for determining 
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view factors in the field of radiation heat transfer. The method described 
herein can be interpreted as though the measured data points are "radiating" 
information to the interpolation point. The farther the data point is away, 
the weaker the "radiation" (lower weight given to the information). In 
principle, the method can easily be extended to any number of independent 

variables, N. 

Step 4. Fit a Polynomial Through the Interpolation Points and Make 

Prediction. 

The final step in the process involves fitting a polynomial through the 

ten interpolation points and then using the polynomial to make a prediction 
of the dependent variable at the target point. The polynomial describes how 
the dependent variable behaves as a function of distance along the prediction 
vector. By trial and error it was found that a forth-order polynomial worked 
well. The polynomial could be used for interpolation or extrapolation 
depending on the location of the target point. There would of course be 

considerably more uncertainty in the prediction f or the case of 
extrapolation. Errors in the ten interpolation points tend to get smoothed by 
the polynomial. 

Reliability Diagnostics of the New Interpolation/Extrapolation Technique 

The method proposed herein provides diagnostics to help assess the 

accuracy of the prediction. The computer program can provide the user with 
averages of the independent variables of the data currently in the database. 
If the independent variables associated with the target point are close to 
the database averages then the user can expect a more reliable result to be 
produced. The coefficient of determination of the polynomial fit through the 
ten interpolation points can be presented to the user to assess the scatter 
in the data. Finally, as shown schematically in Fig. 3.4, the ten 
interpolation points, the polynomial curve, and the prediction can be 
graphically illustrated on the computer screen to show how the dependent 
variable behaves as a function of distance along the prediction vector and 
also indicate the location of the target point along the prediction vector 

relative to the min and max points. If the dependent variable oscillates 

wildly, then unreliable predictions can be expected, particularly in the case 
of extrapolation. In most cases, target points located approximately half way 
between the min and max points produce the best results. 

Testing the New Interpolation/Extrapolation Technique 

The proposed method was tested by making predictions based on a database 
created using a known function so that the reliability of the prediction 
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could be assessed. The form of the function used was: 


D.= lT(rnt + 

1 j=l 


rn# 


x. .+ rn# * .) 

J,i J.i 


(3.3) 


A set of 15 random numbers, rn#, was required - three for each of the 5 
independent variables used to generate a data set. The random numbers were 
held constant during the data set creation so that a consistent set of 
dependent variables were generated. The intent here was to develop an 
unbiased, sophisticated and consistent set of data to provide an objective 
test of the interpolation/extrapolation technique. 

The values of the independent variables used were obtained f rom a 
hypervelocity impact data set. Table 3.1. This data set was selected for two 
reasons. First, it seemed desirable to use a set of actual engineering data 
to provide a realistic test of the technique. Secondly, as is discussed 
below, the technique proposed herein did a poor job of predicting the 
behavior of some of the dependent variables of Table 3.1. Accordingly, it was 
of interest to determine if the nature of the data or the prediction 
technique was at fault for the poor predictions. 

Eq. (3.3) was used with two sets of random numbers to generate two sets 
of consistent data. An analysis was done with each set of data as follows. 
Each record (data point) was temporarily removed from the database, a 
prediction made for the independent variables associated with that record, 
and then the record was returned to the database. This was done for all of 
the 35 records in the database. Thus, all predictions were made with the 
actual data point of interest removed from the database. 

A typical set of results is shown in Fig. 3.5, where actual data values 
are plotted against their corresponding predictions. The dashed line in the 
figure is a linear least squares fit through the data of the form y=mx+b, 
where y is the prediction, x is the actual function value, and m and b are 
parameters to be fit. The coefficient of determination of this fit was 0.937. 
Assuming a functional form of y=x produced a coefficient of determination of 
0.934. Similarly, the other consistent data set produced coefficients of 
determination of 0.961 and 0.959, respectively. Ideally, the prediction, y, 
should exactly equal the actual value, x, which would result in a coefficient 
of determination of unity for the line y=x. 

These results seem to be quite good considering that the dependent 
variables were reasonably complicated functions of fifteen random 
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fable 3.1 Experimental Data from Hypervelocity Impact Tests 


Test 

ID 

Bumper Pr. Wall Proj. 

Thick. Thick. Diam. 

Tb(mm) Tpw(mm) Dp(mm) 

Impact 

Angle 

#(deg) 

Proj. 

VeL 

V(km/s) 

Bump. 

Maj.Ax. 

(mm) 

Bump. 

Min.Ax. 

(mm) 

MLI 
Pen. 
(cm ^2) 

MLI 

Per/Chr 

_(cm_~2).._ 

Pr.Wall 

Maj.Ax. 

(mm) 

Pr.Wall 
Min. Ax. 
(mm) 

221 A 

0.81 

1.60 

6.35 

45.00 

5.52 

15.24 

11.68 

32.26 

51.61 

13.46 

9.65 

227B 

0.81 

1.60 

6.35 

45.00 

7.12 

15.24 

10.92 

64.52 

425.81 

25.40 

12.70 

333 

1.02 

3.18 

4.75 

45.00 

2.88 

10.16 

7.62 

12.90 

32.26 

0.00 

0.00 

334 

1.02 

3.18 

4.75 

45.00 

3.61 

10.16 

7.87 

8.39 

63.23 

0.00 

0.00 

221C 

1.02 

3.18 

4.75 

45.00 

4.57 

11.43 

9.14 

19.35 

51.61 

0.00 

0.00 

22 IB 

1.02 

3.18 

4.75 

45.00 

5.89 

13.72 

10.67 

12.90 

141.94 

0.00 

0.00 

221 A 

1.02 

3.18 

4.75 

45.00 

6.36 

12.19 

10.16 

12.90 

148.39 

0.00 

0.00 

336 

1.02 

3.18 

6.35 

45.00 

4.47 

13.46 

10.67 

64.52 

129.03 

14.99 

6.35 

201B 

1.02 

3.18 

6.35 

45.00 

5.51 

13.46 

10.92 

64.52 

135.48 

13.21 

11.68 

201C 

1.02 

3.18 

6.35 

45.00 

7.21 

13.46 

6.35 

16.13 

161.29 

2.54 

2.54 

203B 

1.02 

3.18 

7.62 

65.00 

3.67 

22.10 

11.94 

22.58 

290.32 

0.00 

0.00 

203A 

1.02 

3.18 

7.62 

65.00 

6.45 

23.88 

13.46 

29.03 

232.26 

0.00 

0.00 

003A 

1.02 

3.18 

7.95 

45.00 

6.51 

19.30 

13.72 

32.26 

129.03 

76.20 

38.10 

338 

1.02 

3.18 

7.95 

45.00 

6.98 

21.34 

14.48 

83.87 

58.06 

25.40 

17.78 

337 

1.02 

3.18 

7.95 

45.00 

7.00 

19.56 

13.21 

64.52 

90.32 

27.94 

12.70 

203F 

1.02 

3.18 

8.89 

65.00 

3.04 

24.89 

12.45 

13.55 

270.97 

0.00 

0.00 

339 

1.02 

3.18 

9.53 

45.00 

6.49 

21.08 

17.53 

129.03 

258.06 

50.80 

38.10 

218B 

1.02 

4.78 

8.89 

45.00 

6.40 

20.32 

15.24 

51.61 

483.87 

18.29 

15.49 

218C 

1.02 

4.78 

8.89 

45.00 

6.76 

21.34 

14.99 

70.97 

270.97 

30.73 

10.16 

230B 

1.60 

3.18 

4.75 

45.00 

3.23 

11.94 

9.14 

3.87 

6.45 

0.00 

0.00 

230A 

1.60 

3.18 

4.75 

45.00 

4.41 

12.19 

9.91 

6.45 

48.39 

0.00 

0.00 

301 

1.60 

3.18 

6.35 

45.00 

2.95 

13.72 

10.92 

4.26 

118.32 

0.00 

0.00 

205A 

1.60 

3.18 

6.35 

45.00 

4.11 

15.49 

12.19 

11.61 

116.13 

5.08 

5.08 

205B 

1.60 

3.18 

6.35 

45.00 

4.59 

16.51 

12.45 

24.52 

335.48 

5.08 

5.08 

205C 

1.60 

3.18 

6.35 

45.00 

5.30 

15.24 

12.70 

16.13 

83.87 

7.62 

7.62 

209B 

1.60 

3.18 

6.35 

65.00 

6.40 

22.10 

13.21 

5.16 

193.55 

0.00 

0.00 

209D 

1.60 

3.18 

6.35 

65.00 

7.40 

19.56 

14.48 

19.35 

206.45 

0.00 

0.00 

207A 

1.60 

3.18 

7.62 

65.00 

5.86 

22.35 

14.99 

11.61 

329.03 

4.06 

4.06 

207C 

1.60 

3.18 

7.62 

65.00 

7.08 

25.91 

16.26 

103.23 

174.19 

0.00 

0.00 

002B 

1.60 

3.18 

7.95 

45.00 

6.39 

20.57 

15.75 

129.03 

77.42 

5.59 

5.59 

211B 

1.60 

3.18 

8.89 

45.00 

5.85 

21.84 

17.27 

77.42 

122.58 

27.94 

12.70 

210B 

1.60 

3.18 

8.89 

65.00 

5.70 

28.70 

16.76 

41.94 

96.77 

3.18 

3.18 

210D 

1.60 

3.18 

8.89 

65.00 

6.80 

35.56 

17.27 

45.16 

212.90 

5.84 

5.84 

303B 

1.60 

4.06 

7.95 

45.00 

4.34 

18.03 

14.48 

32.90 

362.26 

0.00 

0.00 

303 

1.60 

4.06 

7.95 

45.00 

4.59 

18.54 

14.73 

17.03 

166.84 

0.00 

0.00 
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Fig. 3.5 Plot of actual VS predicted values using consistent data. 
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coefficients with five independent variables. 

Applying the New Technique to Hypervelocity Impact Data 

Personnel from the Structures and Dynamics Lab of Marshall Space Flight 
Center (MSFC) provided the author with a set of experimentally obtained 
hypervelocity impact data. Table 3.1. These impact tests were made with the 
multilayer insulation (MLI) placed directly against the pressure wall. The 
bumper plate was placed approximately 100 mm in front of the pressure wall 
plate. For this series of data, the pressure wall was unstressed. As listed 
in Table 3.1 and illustrated in Figure 1.1, the dependent variables measured 
included the major and minor axis dimensions of the bumper hole, the area of 
the hole clean through the MLI called the penetration area, the area of MLI 
outside of the penetration area obviously damaged by the impact called the 
perf orated/charred area, and the major and minor axis dimensions of the 
pressure wall hole. Some comments will now be made on the characteristics of 
these dependent variables. 

The bumper plate hole typically takes the form of a single, well 
defined, relatively smooth, elliptical hole. The greater the impact angle, 
the more elliptical the hole. It is not surprising that the bumper plate hole 
data is the most consistent of all the data given the relatively simple 
nature of the damage. 

The remainder of the dependent variables are much more affected by 
characteristics of the fragmentation/vaporization process of the projectile 
than the bumper hole is. Launch loads typically cause the soft aluminum 
projectile to deform into a variety of nonspherical shapes. This effect, and 
the inevitable presence of a random assortment of microscopic flaws in the 
projectile and bumper, can cause large variations in the nature of the 
particles (from both the projectile and the bumper) that leave the back face 
of the bumper after the bumper-projectile impact. Thus, similar testing 
conditions can produce significantly different damage to the the MLI and the 
pressure wall. 

There is a great deal of inconsistency in the MLI data. In addition to 
the random processes discussed previously, the inconsistency could be partly 
due to the difficulty in visually measuring the areas of damage (penetration 
and perf orated/charred ) because of the rough, irregular shapes of these 
areas. 

Damage to the pressure wall typically consists of a large number of 
craters of various sizes, and possibly some penetrations. The craters and 
penetrations are typically distributed over a relatively large area as can be 
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seen in the photographs of Ref. 3. The data in Table 3.1 gives the dimensions 
of the largest penetration in the pressure wall, which would essentially 
depend on the the largest fragment that results from the bumper-projectile 

impact. As has been discussed, the same test conditions could produce a large 

variation in the size of the largest f ragment and hence the size of the 
penetration. This leads to scatter in the pressure wall data. 

The procedure described previously that was used to test the 

interpolation/extrapolation technique with the consistent data was also used 
with the experimental data of Table 3.1. Each record (data point) was 

temporarily removed from the database, a prediction made for the independent 
variables associated with that data point, and then the data point was 

returned to the database. The predicted versus actual data are shown in Figs. 
3.6-3.11. Also drawn on these figures are solid lines indicating the ideal 
case of ff predicted"="measured" The coefficients of determination associated 
with these predictions are given in Table 3.2. As can be seen from Table 3.2, 
the predictions for the bumper plate are acceptable. The predictions for the 
pressure wall are marginal, although Figs. 3.10 and 3.11 are somewhat 
pessimistic looking since a large number of good predictions were made f or 

data located near the origin of the plots (no penetration case). The 
predictions of MLI damage are poor for penetration area and terrible for 
perforated/charred area. 

Since the proposed interpolation/extrapolation technique produced 
acceptable results for both the consistent test functions of Eq. (3.3) and 

the bumper plate data of Table 3.1, the poor predictions of MLI and pressure 
wall damage are probably due to the scatter in the data produced by such 
effects as the distortion of the projectile during launch and the apparently 

random assortment of microscopic flaws in the projectile and the bumper. The 
proposed interpolation/extrapolation technique appears to be a useful tool 
for engineering design work. 
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Fig. 3.6 Plot of measured VS predicted data for the major axis of the 
bumper plate hole. 
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Fig. 3.7 Plot of measured VS predicted data for the minor axis of the 
bumper plate hole. 
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Fig. 3.8 Plot of measured VS predicted data for the MLI penetration area. 
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Fig, 3.9 Plot of measured VS predicted data for the MLI perforated/charred 
area. 
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Fig. 3.10 Plot of measured VS predicted data for the major axis of the 
pressure wall hole. 
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Fig. 3.11 Plot of measured VS predicted data for the minor axis of the 
pressure wall hole. 
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Table 3.2 Coefficients of Determination for Predictions 



Coefficients of Determination 
(y = prediction, x = measured) 

Data Set 

Line of Form y = mx +b 

Line of Form y=x 

Bumper Major Axis 

0.815 

0.811 

Bumper Minor Axis 

0.774 

0.773 

MLI Penetration Area 

0.322 

0.289 

MLI Perforated/Charred Area 

0.042 

-0.445 

Pressure Wall Major Axis 

0.541 

0.538 

Pressure Wall Minor Axis 

0.575 

0.566 
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4. The Polynomial Function Prediction Technique 

In this section the polynomial function prediction technique is 

described. This method is based on the concepts associated with the finite 

element method (FEM). In FEM, relatively low order polynomials are used to 
interpolate the functions of interest (such as displacements, temperatures, 
and velocities) over a small portion of domain where the function is active 
called an element. The coefficients of the polynomial are derived from known 
values of the function of interest at points called nodes on the boundary of 
the element. For this application, the nodal values of the functions of 

interest (bumper hole size and so forth) were measured experimentally and are 
thus known quantities. This technique involves randomly selecting a 
sufficient number of experimental data (node) points and then determining the 
the coefficients of the polynomial from this data. 

During the course of this project, a more sophisticated polynomial 
interpolation approach was attempted using the isoparametric function mapping 
technique of FEM. This approach in its current state was not found to be 

suitable for engineering trade study purposes. The interested reader can 
consult Appendix 1 for more details on this approach. 

Ideally, the nodes "closest" to the prediction point in impact parameter 
space should be used to evaluate the polynomial coefficients and thus make a 
prediction. However, the set of closest nodes may not form linearly 
independent set of data, making it impossible to solve for the polynomial 
coefficients. Also, one of the closest nodes may have a large experimental 
error which would contaminate the prediction. Accordingly, the computer 
program randomly picks subsets of data from the experimental results database 
file and attempts to make a prediction. If the data is linearly independent 
and a prediction is obtained, then the prediction magnitude and the mean 
distance of the data points from the prediction point are shown on the 
screen. The program is currently designed to seek five predictions. 

In general, the impact parameters will vary greatly in magnitude. In 
hypervelocity impact work, dimensions can be of order 10 and velocities of 
order 10^. This polynomial function approach requires a reasonable scheme for 
determining "distances" between data points in impact parameter space. This 
is accomplished in the program by scaling the impact parameters such that 
their mean value is equal to unity. Of course, the dependent variables, such 
as bumper hole size, need not be scaled. Having scaled the independent 
variables, the usual formula for determining the "distance”, R., between two 
points (experimental data point and the prediction or interpolation point) in 


32 



a multidimensional space can be used: 


- l (x j,i- x j.INT ) ‘ 


The need for scaling the independent variables is evident from considering 
the form of Eq. (4.1). 

The form of the polynomial will now be considered. FEM theory dictates 
that a "complete" polynomial should produce the best results.* 0 Here we have 
four independent variables, x. # . (j = 1 to 4), associated with the i-th 

experimental data point to consider (bumper thickness and so forth). It was 
decided to use Ax^.,. values in the polynomial equation to simplify the 
calculations. The lowest order complete polynomial for this case is: 


C 1 * C 2* Ax l’i * C 3 -ilx 2'i * C 4' ix 3M * C s” Ax 4'i * 

V ix l'i ,4x 2’i * C / 4x ri* 4x 3’i * C s’ Ax 1m‘ Ax 4M * 
C 9 -4X 2 ,.*4X 3 .. * C, 0 *Ax 2 , , -ix: 4 . j 4 C u -4x 3 ,.*ta 4 ,. * 

C 12* Ax l’i" Ax 2'i' Ax 3-i * C 13* 4x l ’ i* Ax 2’ i * Ax 4’ i * 
C 14* 4 Vi‘ 1X 3’i‘ AX 4-i * C 15' iX 2'i' 4x 3M' Ax 4'i * 
C 16* Ax l’ i* Ax 2'i" Ax 3* i* Ax 4’ i 


Sixteen linearly independent data points, D., are required to determine the 
sixteen polynomial coefficients, C.. Eq. (4.2) allows for a linear variation 
in damage along each coordinate axis in the design space. Obviously, 
allowing for a quadratic variation in the damage would provide a much better 
fit to the data. Unfortunately, a "complete" quadratic function with four 
variables would require 81 linearly independent experimental data points with 
the MLI placement and the material types the same for all the data points. 
Currently, experimental data of this nature is not available. 

Coefficient C^ is the prediction of the damage at the point in the 
design space where the prediction is required since this is the value of the 


polynomial when all Ax 


are set equal to zero. 


prediction parameters, such as bumper thickness, 


If one or 
does not 


more of the 
vary in the 


experimental database file then program POLYMETH will sense this and 


automatically take that variable or variables out of Eq. (4.2). If one impact 
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parameter does not vary, only eight polynomial coefficients need be 
determined and thus only eight linearly independent data points are required. 

If there are less than 16 data records in the experimental database 
file, or if the data the data does not span the impact parameter design 
space, then the number of terms in Eq. (4.2) must be reduced if a solution is 
to be found. By pressing the F2 function key, the user can direct the program 
to seek a solution to the following "simpler" equation: 


D i - C 1 * C 3 ,4 Vi + C 4‘ 1X 3'i * S' 4 Vi * 


C 6‘ i! Vi‘ 4x 2'i * C 7* to l'i' AX 3’i * C 8" iX l'i‘ A Vi * 
C 9* AX 2'i* AX 3’i * C 10' 4X 2'i" A Vi * C ll* AX 3'i’ AX 4'i 


(4.3) 


Eq. (4.3) only requires 11 linearly independent data points to obtain a 
solution. However, Eq. (4.3) is an incomplete polynomial which is 
theoretically less desirable than the complete polynomial of Eq. (4.2). If 
solutions can not be found using Eq. (4.3), then the user may press function 
key F3 to request the computer to seek a solution the "simplest" possible 
polynomial: 


D i - C 1 * C 2 * Ax 1'1 


+ C 3* Ax 2-j + C 4* Ax 3’i * C 5* Ax 4’i 


(4.4) 


Eq. (4.4) only requires 5 linearly independent data points for a solution. 

The computer program will repeatedly select random subsets of data from 
the experimental results database file and attempt to find a solution for the 
polynomial coefficients until five solutions have been found. The predictions 
associated with these solutions and the mean "distances" of the data points 
associated with the solutions are printed to the screen. A weighted average 
of the three solutions with the shortest mean distances is also determined 
and printed to the screen. A function of the form of Eq. (3.1) is used to 
determine the weighted average. 

The weighted average should be considered the best value for the 
prediction. If the three predictions with the smallest mean distances are 
consistent then the prediction is probably a good one. 
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5. The Nondimensional Parakcter Prediction Technique. 

In many applications it has been found that empirical functions are best 
represented in terms of nondimensional parameters. Reynolds number is an 
example of a nondimensional parameter that has found widespread use in 
empirical equations of fluid mechanics. Program NONDIMEN uses a series of 
empirical functions based on nondimensional parameters of the form given in 
Ref. 11: 
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(5.4) 


The function coefficients were determined using an optimization routine 
to adjust the values of the coefficients so as to maximize the coefficient of 
determination (R 2 ) of each of the functions. This approach to coefficient 
evaluation is suitable for any form of prediction function - linear or 
nonlinear. The nature of the optimization routine will now be described. 

The magnitudes of the function coefficients can vary by several orders 
of magnitude. To avoid numerical problems it is advisable to work with 
percentage changes in the function coefficients. This approach also provides 
a simple way of controlling the amount of change in the function coefficients 
from one optimization iteration to the next. If the maximum allowable 
percentage change is too large, the optimizer could thrash back and forth 
around the optimum design point without ever converging to it. Alternatively, 
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if the maximum allowable percentage change is too small, then it could take 
an impractical number of iterations to get to the optimum design point, or 

the optimizer could get "stuck" in a local maximum of the coefficient of 

determination function before getting to the global maximum. 

The maximum allowable change in the in the nondimensionalized design 

variable magnitudes is called the "search domain parameter". This is a user 
controlled input parameter. A value of 1.0 (equivalent to a 1007. change) is 

recommended. The optimizer is designed to reduce the magnitude of the search 
domain parameter as the optimization process proceeds. The final value will 
be 1/100 of the initial value. The idea here is to allow large changes in the 
design variables initially, to quickly get into the vicinity of the global 

maximum in the design space, and then use finer steps to precisely locate the 
global maximum. The user is free to change this parameter to attempt to 

improve optimization efficiency. 

The initial values of the function coefficients are set equal to zero. 

Optimal values of the function coefficients could be positive negative or 
zero. 

The method chosen here for search vector selection is based on Powell's 
12 

method .This is a first order method that does not require the calculation 

of the gradient vector. Here, Powell's method was modified as follows. 
Initially, a number of search vectors equal to the number of function 
coefficients are created. The components of these vectors are random numbers 
between -1 and +1. The components of each random search vector are then 

scaled, such that the largest component has a magnitude of unity. These 
vectors are stored as columns of a "search matrix". Next, the coefficient of 
determination is evaluated at the current point in the design space and at 
design points given by +/- the search domain parameter times the first column 
of the search matrix. If either of the + or - design points has a coefficient 
of determination greater than that of the current design point, then the 

design point corresponding to the the highest coefficient of determination 
will become the new design point. Otherwise, the design point does not 

change. The search vector multiplier (+/- search magnitude parameter or zero) 
used with the search vector is stored for later use. This procedure is then 

repeated with the remaining columns of the search matrix. 

A new search vector is created after using all of the search vectors in 
the search matrix. This new vector is created by vectorially adding together 
all of the search vectors times their search vector multipliers. The new 
search vector is a vector sum of previous successful search vectors since 
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unsuccessful search vectors have search multipliers of zero. Thus, the new 
search vector represents (stores) the trend of the optimization process. The 
new search vector is scaled such that the magnitude of it’s largest component 
is unity and then is used to replace the first column of the search matrix. 
The procedure is repeated, a new search vector is determined, and then used 
to replace the second column of the search matrix, and so forth until only 
the last column of the search matrix remains untouched. Then an entirely new 
search matrix is created using the random number generator, and the process 
continues. 

If at any time in the iterative process, a new search vector has a 
magnitude of zero (implying all current search directions are not 
beneficial), then a new random search matrix is created immediately. The 
random number generator uses a seed based on the number of seconds from 
midnight on the computer’s clock. Each successive rim of the optimizer will 
use a different set of search vectors. Currently, the program runs the 
optimizer three times (each time using different sets of random search 
vectors) to help ensure that the global maximum of the coefficient of 
determination has been located in the design space. 

The number of search search matrices generated is governed by a user 
input parameter called the "iteration parameter". The number of random search 
matrices generated is equal to the number of design variables times the 
iteration parameter. The recommended value for the iteration parameter is 20. 

As can be seen from the test runs of Table 5.1, the optimizer produced 
very consistent coefficients of determination for all four prediction 
equations (Eqs. 5. 1-5.4). It was noted that virtually identical coefficients 
of determination could be produced by prediction functions having very 
different coefficient magnitudes as is illustrated in Table 5.2. This is a 
typical characteristic of nonlinear equations. 

After the prediction function coefficients have been determined and 
displayed on the screen, the user will be prompted for the impact parameters 
(such as bumper thickness) associated with the desired predictions. Multiple 
predictions can be made from the same set of prediction coefficients. 
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Table 5.1 Prediction Function Coefficients of Determination for Several Runs of the Optimizer 
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Table 5.2 MLI Hole Diameter Prediction Coefficients for Several Optimizer Runs 
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6. A Comparison of the Accuracy of the Prediction Techniques 

The accuracies of the three prediction techniques discussed herein were 
compared with respect to a common impact data set. Table 6.1. This is the 
same data set that was recently tested for insulation damage in the Sunspot 
Thermal Vacuum Chamber. This data is also provided on the computer disks as 
experimental database file MLI.DAT. These specimens had the MLI mounted next 
to the bumper during impact testing. Ref. 11 contains more general details on 
the impact testing. 

The accuracy of each prediction technique was tested by first removing a 
data record from the experimental database file, and then using the remaining 
data to make a prediction for the impact damage associated with the impact 
parameters of the removed data record. This was repeated for all of the 19 
data records of Table 6.1. The results of this accuracy check are shown in 
Tables 6.2 to 6.4 for the three prediction techniques. To compare the 
accuracies of the three prediction techniques, average percentage differences 
were calculated f or each of the f our prediction f unctions. These are 
summarized in Table 6.5. Here, average percentage difference is the average 
magnitude of the difference between the predicted and measured values divided 
by the average magnitude of the measured values, times 100. Thus, relatively 
high average percentage differences indicate that the prediction function did 
a poor job of predicting the damage. 

The following observations can be made about Table 6.5: 

1. The poorest predictions by far were made for the pressure wall hole 
diameter. 

2. The best predictions were made for the minimum bumper hole diameter. 

3. The inverse R and nondimensional functions did an acceptable job for 
engineering trade study purposes (average percentage differences < 207.) for 
predicting the bumper hole size and the MLI hole diameter. 

4. The nondimensional function technique did the best job overall of 
predicting impact damage. 

The nondimensional function approach did the best job of predicting the 
data of Table 6.1. However, different data sets could produce significantly 
different results. The nondimensional function approach may not work as well 
if the prediction parameters (such as impact velocity) cover a greater range 
in the database. Also, the inverse R method has the advantage of being able 
to easily incorporate additional impact parameters. The other two prediction 
techniques are not as flexible. For instance, the inverse R method would be 
the method of choice for the case where different materials are used for the 
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1032 1.60 3.18 4.75 45 4.75 12.12 10.19 25.40 0.11 0.00 0.00 

1024 1.60 3.18 4.75 0 3.79 10.97 9.53 17.78 0.15 0.00 0.00 

1031 1.60 3.18 6.35 0 6.98 13.69 12.73 38.10 1.05 0.00 0.00 



Table 6.3 Typical Prediction Results for Polynomial Impact Damage Function 
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Averages! 0.514 0.511 I 0.561 0.563 1 1.637 1.752 I 0.113 0.272 

Ave % Dif 10.04 16.23 27.22 158.27 



Table 6.2 Typical Prediction Results lor Inverse R Impact Damage Function 
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Averages J 0.514 0.515 I 0.561 0.550 I 1.637 1.616 1 0.113 0.113 

Ave % Dif 8.87 15.15 15.99 92.69 


Table 6.4 Typical Prediction Results for Nondimensional Impact Damage Function 
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Averages I 0.514 0.518 I 0.561 0.562 I 1.637 1.587 I 0.113 0.131 

Avc % Dif 5.42 6.21 10.88 66.65 
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bumper in the same experimental results database file. 
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7. Conclusions and Recommendations 

As a result of this study the following conclusions were reached: 

■ There is a large amount of scatter in the hypervelocity impact damage 
data. It is doubtful that very high prediction accuracies can be obtained 
regardless of the prediction technique used. 

■ There is not a great deal of data available for any given set of impact 
conditions (such as the case with MLI against the bumper). Lack of data 
prevents higher order prediction functions from being used. 

■ The inverse R method is the most flexible prediction scheme. Any number of 
impact parameters and any size of database can be treated. 

■ The nondimensional parameter functions seem to do the best job of 
predicting impact damage over a relatively restricted range of impact 
parameters. 

Based on this study the following recommendations are made: 

■ If possible, all three prediction techniques should be evaluated to 
determine the best possible prediction technique for a given data set. 

■ The nondimensional parameter scheme should be used to make impact 
predictions from data sets for which the impact parameters have a relatively 
small range. 

• The inverse R prediction technique should be used in applications where 
there are a large number of impact parameters (different bumper materials in 
a single database file for instance) or where the impact parameters vary over 
a wide range. 

■ Numerical simulation results or approximate analytical results for high 
velocity (10 - 15 km/sec) should be placed in the "experimental" results 

database file so that realistic predictions for on orbit impacts can be made 
with the software. 
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Appendix 1 A Prediction Technique Based on an Extension of the 
Isoparametric Formulation of the Finite Element Method - A Master’s Thesis by 
Mr. P. Wang. 
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CHAPTER ONE 


INTRODUCTION TO MULI VARIABLE ANALYSIS 


1.1 Multivariable Analysis 

Multivariable analysis is concerned with data that is a function 
of several Independent variables, which is very important and common 
not only in engineering analysis but also in other fields [11. Many 
multivariable analyses are derived from the viewpoint of statistical 
theory and have been already well developed such that they can help 
engineers make empirical predictions from their results. Regression 
analysis [2] is a technique that is commonly used to analyze 
experimental data in various areas of research. 

Interpolation schemes can provide powerful tools for determining 
the relationship between dependent and independent variables. For 
example, interpolation of scattered 2-dimensional and 3-dimensional 
data using the Shepherd method, has been an important subject of CAGD 
(Computer Aided Geometric Design) recently [31. Many different 
Interpolation schemes have been devised for various types of scattered 
data [4]. Most of these methods such as data point triangulation and 
B-spline interpolation, are in essence a type of surface fitting [51. 
Presently, attention is focused on how to reduce geometric 
discontinuities, smooth the data error and apply these procedures in 
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computer graphics [6]. Tadeusz Liszka has proposed a local 
interpolation method by using a Taylor expansion of the unknown 
function to reduce geometric discontinuity for those schemes that fit 
scattered data [7]. 

Due to the rapid development of computers in recent years, the 
finite element method (FEM) has become an important tool for the 
solution of engineering problems. FEM is a numerical approximation 
procedure based on interpolating the variables of interest over finite 
parts of the continuum called elements. The isoparametric formulation 
is one of the more important implementations of FEM, where the element 
coordinates and element displacements are both interpolated using the 
same shape functions that are defined in a natural coordinate system 
[ 8 ]. 

However, currently most of the practical applications of CAGD and 
FEM are focused on two dimensional or three dimensional problems, and 
little attention is paid on the cases with more than three dimensions 
or independent variables. 

1 . 2 Purpose of Study 

The original motivation of this thesis was to derive a 
formulation from the isoparametric concept of FEM to predict the damage 
to spacecraft in low earth orbit from space debris traveling at 
hyperve loci ties. Because a number of physical and mechanical 
properties of the debris and spacecraft are expected to be related to 
the damage, a multivariable analysis is thus required to make 
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predictions from experimental data. Although Bouma and Burkitt and 
several others have worked on this problem since 1963, the methods used 
were typically based on statistical model theory [9]. The author was 
interested in extending the isoparametric concept of FEM to construct a 
new model for use in multivariable analysis. There are many different 
areas of research, such as geostatistics in the field of geology [10] 
and biostatistics in the field of biology [11], which require 
multivariable analyses. The FEM based model developed here was 
designed to be applicable to multivariable analysis problems in 
general . 

The following chapter reviews the basic concepts of the 
isoparametric implementation of FEM and then extends the concepts to 
problems with an arbitrary number of degrees of freedom (DOF]. Chapter 
3 discusses the influence of element distortion on calculated results. 
In Chapter 4, two sets of actual experimental data are used to test 
the proposed model. Here, the results of the new model are compared 
with those from linear regression. Conclusions and recommendations are 
presented In Chapter 5. Listings of the computer programs developed 
for this investigation are given as an Appendix. 
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CHAPTER TTIO 


MODELING USING THE ISOPARAMETRIC CONCEPT 


2.1 Generalization of the Isoparametric Concept 

The finite element method Is basically a discretization process 
to partition a complicated structure or system Into a finite number of 
small parts having simple geometric shapes [12]. These small parts are 
called elements. A group of elements modeling a continuum is called a 
mesh. Points In the continuum at the corners (and sometimes along the 
edges) of the elements are called nodes or nodal points. The FEM 
explicitly determines values for the dependent variables at the nodes. 
Simple functions are chosen to approximate a physical quantity over 
each finite element. Such assumed functions are called interpolation 
functions or shape functions, which are functions with unit value at 
one nodal point and zero value at all other nodal points. Through the 
use of shape functions, a relationship can be established between the 
coordinates of every point Inside an element and the element nodal 
coordinates (called degrees of freedom, DOF). The principal idea of 
the isoparametric formulation consists in using these same shape 
functions to Interpolate the physical quantity of interest over the 
element [13]. Thus, a similar relationship can be established between 
the physical quantities at every point inside an element, and the 
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element nodal physical quantities. To implement the isoparametric 
formulation, an orthogonal natural coordinate system is introduced such 
that elements described in the physical coordinate system can be mapped 
into an element in the natural coordinate system, where each coordinate 
axis varies from -1 to 1. In order to illustrate this point, a 1-DOF 
case is derived as follows. 

Consider a bar element with two nodes which lies along the 
X-coordinate axis, as shown in Fig. 2. 1. Because we want to have the 
whole element mapped from the physical coordinate X to the natural 
coordinates, say £, where -ls£sl, the following correspondence 
(boundary conditions) must occur : 

1. when £ - -1, X = X i ■, 

2. when £ = 1, X * X 2 

where X i and X are the nodal X coordinates. The shape functions here 
must be linear in form since there is 1-DOF and only two boundary 
conditions. Thus, a suitable equation for writing the physical 
coordinate as a function of the natural coordinate Is : 

X - -5- (1 - £)X i + -|- (1 + £)X z (2.1) 

or 

2 

x - £ N X (2.2) 

1 = 1 

where (1 - £) and N 2 = (1 + £) are the shape functions. 
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Figure 2.1 A 1-DOF bar element with 2 nodes shown with respect to 
the physical coordinate system (X) and the natural 
coordinate system (£)• 


65 


Equation 2.1 or 2.2 establishes our desired mapping relationship. 

Now, if we are given some nodal physical quantity, say displacements 

and LJ , for these two nodes, then based on the principal of 
2 

isoparametric formulation, the displacements U at any position in the 

bar can be determined by the same shape functions N^and That is, 

2 

U = V N U (2.3) 

u \ 1 
1=1 

A similar derivation car be made for the 2-DOF and 3-DOF cases except 
that the shape functions will have different forms. More discussion 
on this topic follows in next section. 

After understanding the basic concept of isoparametric 
formulation, we begin to extend this concept by generalizing the 
geometric coordinates to be any physical coordinates such as, say, mass 
velocity or concentration. 

2.2 Shape Functions for N-dimensional Analyses 

As it was introduced in Section 2. 1, shape functions play a 
paramount role in relating some physical quantity within an element to 
the element nodal values. The form of the shape functions depends on 
the number of DOF of the problem and the number of nodes in the element 
[14]. For the case of a 1-DOF bar element, each node just had one DOF, 
while for the 2-DOF case such as a quadrilateral (four sided) element, 
each node will correspond to 2-DOF., because we need at least 2 
coordinates to describe the "location" of each node in the 2-D plane. 

In this thesis, we will focus on the linear Interpolation in each 
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DOF direction over an n-DOF element. The procedure that was used for 
constructing the shape functions for the 1-DOF element will be 
generalized for use in 2-DOF, 3-DOF and n-DOF elements in the following 
discussion. 

Starting from the 2-DOF element, we need 4 nodes to make an 
element such that a linear interpolation can work along the 2-DOF 
directions associated with the element. This element is called 
bilinear element, as shown in Figure 2.2-a, and the corresponding four 
shape functions in terms of the natural coordinates can be found 

by using an approach similar to that which led to equation (2.1) and 
written as 


N = 
1 

i 

4 

(1-0 

(l+7») 

(2.4) 

N = 
2 

1 

4 

(1-0 

tl-n) 

(2.5) 

N = 

3 

1 

4 

(l+O 

(1-TJ) 

(2.6) 

N = 

4 

1 

4 

(l+o 

(l+7») 

(2.7) 

:i + 


(1 + T 

^Tf) (1=1.. 4) 

(2.8) 


where £ 7 ) are defined to vary between -1 and 1 [15]. In equation (2.8), 
the sign of the £ and tj terms are determined by the coordinates S^and 
7 )^, which can be both i 1 as long as the node number Is assigned. For 
example, for 1=2, € = -1 and f)= -1, so, N = 1/4(1 - £) ( 1 - • 

Based on the bilinear case, it is not difficult to derive the 
general form of the shape functions for a 3-D element with 8 nodes, 
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called a trillnear element (Figure 2.2-b). The general form Is simply 
written as 

N » -S-(1 + £.0(1 + + CO (1=1. .8) (2.9) 

i o i * * 

where £ , V and C are equal to ± 1 when i Is assigned. 

M 1 l 

In practical engineering problems, the application of FEM Is 
limited to 2-DOF or 3-DOF elements because the physical quantities are 
typically assumed to be the function of spatial coordinates. If we are 
given a problem based on physical parameters instead of geometric 
coordinates, and the DOF of each node is more than 3, then a more 
generalized shape function is required to correlate these parameters to 
some physical quantity. It is not difficult to extend shape functions 
for 2-DOF and 3-DOF cases to apply to arbitrary n-DOF cases : 

N = n (-^-) n (l ± X ) (2. 10) 

112 n 

n - number of DOF of each node 

where each of corresponds to one of 2 n "corners" of an n-DOF space, 
and X are defined to vary between -1 and 1 as before. 

n 

To this point, we have established an n-dimensional interpolation 
model by n-dimensional shape functions. Next, we want to apply the 
generalized isoparametric formulation to make multivariable 
interpolation for some practical problems. In order to illustrate the 
application of the shape functions for this objective, an example 
involving damage to spacecraft caused by space debris traveling at 
hypervelocities will be considered. 
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In low Earth orbit (LEO), there exists a large quantity of 
orbital debris that has been generated by man’s activity in space over 
the last three decades [16]. The debris varies greatly in size — from 
essentially intact upper stages of rockets to small particles produced 
by explosions on orbit. The most dangerous particles for active 
spacecraft in LEO are of characterist ic dimension ranging from about 
0.5 mm to 2 cm because there are vast numbers of these particles, they 
have a high energy content, and they are too small to track by radar or 
other means. These particles are traveling at orbital hypervelocit ies 
(10-20 km/sec) and thus can inflict a significant impact damage to 
spacecraft. To reduce the impact damage to a minimum, Whipple [17] 
proposed that; a protective device for a spacecraft, called a bumper, 
which consists of a thin aluminum outer shell or plate placed some 
distance from the main hull (pressure wall) of the spacecraft, Figure 
2.3. The function of the bumper Is to break-up or vaporize the debris 
particle. The pressure wall is then impacted by many tiny particles 
rather than a single large one. As a result, the main hull of the 
spacecraft sustains little damage. 

In order to ’lnderstand the damage that such particles will 
inflict upon a spaceship, an experimental approach was developed at 
Marshall Space Flight Center to study the damages produced in simulated 
spacecraft targets by projectiles fired at hypervelocities [18]. A 
large amount of experimental data has been collected over a wide range 
of impact conditions. In this thesis, the author is interested in 
developing a new technique for predicting spacecraft damage from the 
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characteristic diameter of debris 



0=impact angle of debris 



pressure wall 


Figure 2.3 The bumper designed for shielding the pressure wall of the 
spacecraft from hypervelocity impact of orbital debris. 
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experimental results. Many attempts have been made by others to use 
conventional techniques to fit functions to the data with mixed 
success. 

For illustration purposes, a simplified 2-DOF version of the 
impact problem will be considered first (Figure 2.4). Here it is 
assumed that the damage D^is only a function of the velocity of debris, 
v , and the thickness of bumper, t ^ DCv^t^). could for 

instance be the size of the hole produced in the pressure wall by the 
impact. The experimental data would consist of a number of data points 
of the form : 


D = D 

( 

V , 

t ) 

l 

1 

i 

D = D 

( 

V , 

t ) 

2 

2 

2 

d 3 = d 

( 

V 

V 


A prediction of the damage D at some v = v and t * t is 

required. By the concept developed in Section 2.1, we can treat these 

independent parameters, v^and t f as physical coordinates. Before we 
apply the shape functions, a procedure to select the most appropriate 
set of 4 "nodes" is required to make a prediction. These four nodes 
are considered to make up a finite element. This node selection 

procedure will be discussed in next section. Assuming the appropriate 

nodes have been chosen, the equations relating (v,t) in the design 
space to the nodal values are : 
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( 2 . 11 ) 


4 

v = 7 N v 
u i i 
\ =1 

4 

t = 7 N t (2. 12) 

u i l 
i = i 

• • 

where v , are known nodal values. By setting v = v and t = t , 
there are two nonlinear simultaneous equations with two variables, £ 
and 7 j ( to be solved because N^ are functions of natural coordinates. 
There are no direct methods available for solving these nonlinear 
simultaneous equations. Newton’s method, which is based on truncating 
the Taylor series to only linear terms can effectively transform a 
nonlinear system of equations a linear system [19]. In general, an 
iterative approach must be used, which requires an initial guess. More 
details on the influence of the initial guess on uniqueness of the 
isoparametric mapping process will be given in next chapter. A 

subroutine called ’ nonlinsol’ was written to use Newton’s method to 

• • * 

find the roots (appendix A-2). The roots, say £ and , like v and 
t, are proposed to be related to D . Thus, by using the same shape 
functions and substituting and tj* into the equation results in 

D* = £ N D (2.13) 

i=l 

Thus a prediction for D is obtained. 

The same procedure can be followed for the 3-DOF case. If one 
more parameter is added, say 0, the impact angle, then we have D = 
D(v,t,0), and after choosing a set of 8 "nodes", we have 
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(2.14) 


8 

v = y n v 

u i i 
1 = 1 

8 

t ■ ? H t (2. 15) 

u i i 
i=i 

8 

0 = V N 0 (2.16) 

u l i 
i = i 

By solving 3 nonlinear simultaneous equations, 3 roots will be 
obtained to get D like the 2-DOF case. 

Therefore, when we have an n-DOF case, that is, we have n 
independent variables like v,t,0, the same logic is repeated, but the 
number of DOF increases. Of course, the correspondence of ’ s and the 
2 n corners of a linear element in an n-DOF natural coordinate space is 
beyond the geometric imagination of the human mind. A Pascal 
subroutine named ’ shpsign* (see appendix A-l) was written to generate 
the shape functions for an n-dimension element. 

2.3 Strategy of Choosing Element Nodal Values 

In the example of section 2.2, we outlined the basic procedure of 
how to perform multivariable interpolation using the isoparametric 
concept with the generalized shape functions. Here we discuss the 
problem of how to select the most appropriate nodes for a prediction. 

Consider the 2-D example of D = D(v,t) of the previous section. 

• * 

Our goal is to predict D , which is associated with known parameters v 

and t , based on the set of known Referring to Figure 2.4, we can 

* * 

see all the points with the coordinates of (v .t^), including (v ,t ), 
are scattered on the v-t plane. Based on the isoparametric concept of 
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FEM, we need 4 nodes to make a bilinear element, such that we can 
interpolate D* at (v*,t*) based on the 4 known nodal However, 

there exists many possibilities of combining 4 sets of nodes to form a 
bilinear element among the scattered data. It Is reasonable to choose 

the 4 nodes that are "closest" to (v ,t ). Therefore, we have to 

• * 

determine all the "distances" from each (v .tj) to (v ,t ), and sort 
them by the order from the closest to farthest. Typically, the 
physical coordinates will vary greatly in magnitude. For the example 
considered here, the velocities are of order 10 and thickness of order 
1. Thus some form of scaling is required before "distances" from point 
to point in the design space can be determined. The mean value of each 
coordinate was used as scaling factor here. Thus the distance between 
(v .tj) and (v , t ) is : 



number of sets of experimental data. 

Obviously, it is easy to extend this formula to the n-D case — 
the number of the nodes required to form a linear element will become 
2 n . Theoretically, we expect a reasonable approximation made by using 
the closest 2 n nodes. However, using the closest 2 n nodes will not 
necessarily produce the best predictions because of the possible 
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Influence of element distortions. Therefore, on some occasions, we need 


to change one or two, or even all the nodes to reduce the geometric 
distortion of the element. This requires that different elements 
formed by different sets of nodes be tested. More attention will be 
paid on this problem in next chapter. 
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CHAPTER THREE 


THE INFLUENCE OF ELEMENT DISTORTION ON CALCULATED RESULTS 


3,1 Element Distortion Caused by Improper Node Numbering 

The influence of element distortion has been an important subject 
in FEM, because distorted elements may produce poor results [20], As 
stated in Section 2.1* the isoparametric formulation is a one-to-one 
mapping from a set of global cartesian coordinates to a set of local 
(natural) coordinates. Highly distorted elements corrupt the mapping 
process thereby producing unreasonable results [21]. There are two 
possible sources of element distortion : 

(1) distortion caused by Improper elemental node numbering. 

(’2) distortion caused by the geometric Irregularities of the elements. 

The latter source of distortion Is discussed in the next section. 
The first type of distortion may be cured by proper node renumbering 
such that a non-twisted element can be obtained, as shown schematically 
for a 2-D case In Figure 3-1. In order to explain this, a 2-D case is 
considered. Figure 3.2 shows four nodes that are numbered by 1-2-4-3 
in a counterclockwise (CCW) sense such that a bow-tie element Is 
obtained. Based on the concept of isoparametric mapping, this twisted 
element defined In the physical coordinate system (X, Y) is mapped onto 
an element defined in the natural coordinate system Hence, all 
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Physical Coordinate System 


Figure 3 


2 An example of improperly numbered element results in a 
nonunique mapping. 
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the points, including the four corners. Inside the twisted element 
should be uniquely mapped onto the points of a square element in £-tj 
plane. It is easy to investigate the uniqueness of the mapping by 
checking if the corners of the twisted element uniquely correspond to 
those of the square element by the node number. For example, the 
coordinates of node 1 in X-Y plane is (2,4), which is supposed to 
correspond to the coordinates of (-1,1) of the node 1 in £-7) plane, and 
so on. Now, consider corner 3, whose coordinates In X-Y is (5,3) and 
the corresponding coordinates of node 3 In the £-tj plane should be 
(1,-1). By recalling the equations (2.11) and (2.12), two simultaneous 
equations are obtained : 

4 

5 = 7 N X 

u i i 
i = i 

4 

3 = E N Y 

" i l 
i = i 

where i = 1..4 are given by equations (2.4) to (2.7), and X t and Y 
are the X and Y coordinates of the ith node, respectively. After 
solving these two nonlinear equations with the initial guess of (?,tj) = 
(0,0), a set of roots is obtained as £ = 7.0, 7) = 0.08, which Is 
obviously different from the required values of (1,-1). On the other 
hand, if the four nodes are numbered in a cyclic way (1-2-3-4, CCW), 
then a nontwisted element is obtained (Figure 3.3). Here, node 4 has 
(x,y) coordinates (5,3) and should have (£,t») coordinates (1,1). 

Again, by substituting these corners into equations 2.11 and 2.12 and 
solving the nonlinear equations with the initial guess of (0,0), we do 
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a unique mapping. 
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get the root of (1,1) as was expected. The other three corners of this 
element are also mapped to appropriate (£,t?) values. Thus, care must 
be taken to avoid element distortion caused by Improper node numbering. 

A systematic algorithm for determining the proper node numbering 
for a given element will now be discussed. Suppose a set of four nodes 
are given to form an element In X-Y plane (Figure 3.4-a), which must be 
mapped onto an element in plane, whose four corners are already 
cyclically numbered 1-2-3-4, (Figure 3.4-d). First, these four nodes 
are partitioned Into two groups by sorting their X-coordinates from the 
smallest to the largest, Such that the first group contains two nodes 
with the smallest X-coordinates and the other group contains the 
remaining two nodes with the larger X-coordinates. In the first group, 
the node with the least X-coordinate is numbered 1, and the other one 
is numbered 2 ; in the 'Second group, the node with smaller X-coordinate 
is numbered 3 and the other one Is numbered 4 (Figure 3.4-b). Next, 
each group is separately sorted according to the Y coordinates of its 
two nodes. In the first group, the node with the largest Y-coordinate 
is numbered 1 and the other one is thus 2. Similarly, in the second 
group the node with largest Y-coordinate is numbered 4 and the other 
one is numbered 3 (Figure 3.4-c). After these two sortings, a set of 
cyclically-numbered nodes (1-2-3-4) is obtained (Figure 3.4-d). We can 
extend this sorting scheme to a 3-D case, where there will be 8 nodes 
to be numbered to make a least distorted hexahedral isoparametric 
element, Figure 3.5. As was done for the 2-D case, first the nodes are 
sorted with respect to their X components. Then, we partition these 
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Figure 3.5 A 3-D undistorted element is formed by cyclic numbering. 
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nodes into two 4-node groups by the temporary order 1-4 and 5-8. Again 
we treat each group as a 2-D case, and obtain a second temporary order 
for each group by sorting their Y— components. At last, a third sorting 
for Z-components is conducted such that (node l)z > (node 2)z, (node 
4)z < (node 3)z for the group of 1-4, and (node 5)z > (node 6)z, (node 
8)z < (node 7)z for the group of 5-8, so a normal element is obtained 
by the last order. 

The same node sorting technique can be applied to problems with 
more than three dimensions. A pascal subroutine named "nodeswap" is 
written for this purpose (see Appendix 3). 

3.2 Element Distortion Caused by Geometric Irregularities 

In last section, minimizing element distortion by selecting an 
appropriate cyclic node numbering scheme was discussed. However, there 
are often geometric irregularities which cannot be removed by proper 
node numbering. Figure 3.6. These geometric defects may also cause a 
nonunique isoparametric mapping to occur. 

In the FEM isoparametric formulation, the shape functions are 
always used to correlate the element nodal values to the values within 
the element. However, it is possible for our purpose here that the 
points to be predicted are outside the elements, so that the use of 
extrapolation is required. In this case, element distortion could make 
the errors that are inherently associated with extrapolation worse. 

Due to these two considerations, different elements, like 
different meshes in FEM, are tried to find the elements that are the 
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least distorted and that contain the point in space where the 
prediction is required. Thus, in addition to the element formed by the 
closest 2 n nodes, we can try other possible elements formed by 
arbitrarily taking 2 n nodes from the first closest 2 n + 1 nodes to make 
distinct combinations, from which 2 n + 1 different elements will be 
generated. Similarly, elements based on the combinations of the 
closest 2 n + 2 ,etc. nodes also can be further tested. 

Next, how the geometric irregularities influence the uniqueness 
of the mapping will be considered. The mapping from physical 
coordinate space to natural coordinate space requires the solution of 
nonlinear simultaneous equations. In general, nonlinear equations can 
have more than one set of real roots, which can be found by starting 
from different initial guesses when Newton’s method is used. If the 
element has just a small amount of distortion, then one set of roots 
(£,tj) must uniquely exist inside or on the element borders defined by £ 
- ±1 and 7 } = ±1. However, it was found that some highly distorted 
elements still could give a unique mapping for some regions within the 
elements. If the point is located in the vicinity of the less 
distorted part of the element, then the mapping of this point could be 
unique. In order to explain this, an example will be considered in 
Section 3.2. 1. 

If the point to be predicted is outside all possible elements, 
then the natural coordinate roots must be outside the elements as well. 
If the linear interpolation functions are assumed to be valid outside 
the elements when the point of interest is sufficiently close to the 
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elements, then the same concept Just discussed for treating a point 
inside a distorted element is extended to this case. That is, the 
uniqueness of mapping (or the roots of £»7)) depends on if the point of 
interest is located in the vicinity of the undistorted part of the 
element. An example in Section 3.2.2. will be used to explain this 
application. 

3.2.1 Example 1 

In this section, the 2-D "debris impact" example, which was 
introduced in Chapter 2, will again be considered. As shown in Table 
3. 1, there are five nodes defined in the velocity-thickness (v-t) 
coordinate plane, where v refers to the velocity of debris and t is the 
thickness of bumper. Each of the nodes has an associated impact 
damage D, arbitrarily defined for illustration purposes by the 
function: 

D(v,t) = v 2 t + v + t (3.1) 

A prediction of the damage will be made at (v,t) * (3,3) using the 
proposed method with the data of Table 3.1. Prediction will be made 
using five different elements and the results compared with the "exact" 
answer given by equation 3.1. 

Now, we begin the analysis with finding the "distances" from each 
node to the node (3,3) based on the scaling scheme Introduced in 
chapter 2, and then sorting them by the order from the closest to the 
farthest, which is listed in Table 3.2. Next, we arbitrarily select 4 
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Table 3. 



Table 3.1 Data list of example 1 



2 Sorting the data of Table 3.1 by scaled distances 
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sets of nodes from these 5 nodes to make 5 different elements, of which 


the corners are cyclically numbered and then mapped to a natural 
coordinate plane. Each of these elements is separately discussed as 
follows: 

Element 1 is formed by the nodes 2, 3, 4, and 5 in Table 3.2. 
After cyclically numbering the nodes, the shape of the element is shown 
in Figure 3.7-a. Obviously, the node (3, 3) , marked with is inside 

the element. By solving the nonlinear simultaneous system for mapping 
this node to the natural coordinate the unique root (£,tj) = (0.65,0.29) 
is obtained whenever initial guesses between -1 and 1 are provided to 
the equation solver. If the initial guess is exactly one of the 
corners of the element in natural coordinate plane, except the corner 
(1,-1) which results in the roots outside the range of -1 and 1, the 
roots obtained are also (0.65,0.29). However, we are Just concerned 
with the uniqueness of the roots inside the element, so the roots 
outside the element will be ignored. At last, we substitute the roots 
of (0.65,0.29) to equation (2.13) to determine the damage. 

Element 2 is formed by the nodes 1, 3, 4, and 5 in Table 3.2. 

The same approach is repeated as in element 1, but the node (3,3) is 
outside the element as shown in Figure 3.7-b. In the process of 
solving nonlinear system, no matter what the initial guesses are 
between -1 and 1 or the corners, the root we obtain are consistently 
(3,2). By substituting this root to equation (2.13), a damage is 
predicted. 

Element 3 is formed by the nodes 1, 2, 4, and 5 in table 3.2. As 
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shown in Figure 3.7-c, the node (3,3) is still outside the element, and 
the nodes 1, 2, and 5 are co- linear. No convergent roots can be 
obtained whatever initial guesses are chosen because of the serious 
geometric irregularity of the element. Thus, the damage cannot be 
correctly calculated using this element. 

Element 4 is formed by the nodes 1, 2, 3, and 5 in table 3.2. As 
shown in Figure 3.7-d, this is also a triangular element with the nodes 
1, 2, and 5 collinear, but here the node (3,3) is inside the element. 
Except the initial guess by the corner of (-1,1), which causes a 
singular matrix in Newton's method and fails to solve the roots, the 
other guesses uniquely generate (0.38,0.42). Based on this roots, the 
damage can be calculated in spite of the geometric irregularity of the 
element . 

Element 5 is formed by the nodes 1, 2, 3, and 4 in table 3.2. As 
shown in Figure 3.7-e, it is also an irregular element that contains 
the node (3,3). The mapping is not unique because the roots are (0.48 
, 0.24) when the initial guesses are taken from in the area specified 
by 0 s 1 and -1 s i| s 1, but the roots are (-0.7, 1.8) when the 
initial guesses are chosen from the area specified by -1 s ? < 0 and -1 
s 7 ) s 1. Thus damage predictions cannot be made using this element. 

Damage prediction based on five elements discussed previously are 
given in Table 3.3. The prediction made by the first element should be 
most convincing because the element is regular and the node predicted 
is inside the element. The prediction value of 34.8 agrees quite well 
with the exact value of 33, considering that the nodal data values by 
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element 

* 

roots 

* * 

roots 

• 

damage 

■ • 

damage 

2345 

0.65,0.29 

0.65,0.29 

3.48E+01 

3.48E+01 

1345 

3.0 ,2.0 

3.0 ,2.0 

2. 14E+01 

2. 14E+01 

1245 

— 

— 

— 

— 

1235 

0.38,0.42 

— 

3.38E+01 

— 

1234 

0.48,0.24 

- 0.7 ,1.75 

3.41E+01 

2. 96E+01 


note: 

root* : root obtained by initial guess of (0,0) 

root**: root (£,t>) obtained by initial guess of closest corner 
• * 
damage : damage based on root 

• * ** 
damage : damage based on root 

exact damage : 33 

Table 3.3 Comparison of the predicted damages of example 1, based on 
different elements and initial guesses 
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sin order of magnitudes (4 to 43). A discussion on the results of 
element 2 and element 3 will be given after the second example is 
considered. Here, two observations based on the results of element 1, 
element 4 and element 5 are given. 

(1) Elements containing the predicted node are not necessarily 
associated with accurate predicted results : 

Based on providing the equation solver with the initial guess of (0,0), 
it is easy to check that elements 1, 4, and 5 contain the predicted 
node (3,3). However, not all the three elements guarantee a unique 
mapping, which can be seen when other initial guesses are used. 

Element 4 and 5 fail to give reasonable predictions because the 
uniqueness of mapping does not exist for them. Hence, an accurate 
prediction depends on to the uniqueness of mapping if the predicted 
point is known to be inside the element. 

(2) Making predictions based on different elements: 

It is apparent that element 1 can make a better prediction due to its 
relatively undistorted shape which insures the uniqueness of the 
natural coordinate roots inside the element. However, it is difficult 
to determine the geometric irregularity of elements by trying all 
possible initial guesses to test the uniqueness of the mapping. Also, 
uniqueness mapping could occur in highly distorted elements if the 
point to be predicted is located in the less distorted parts of the 
elements. The following approach was used to cope with this problem. 
First, the corner closest to the point to be predicted is located for 
each element in the physical coordinate system. When an isoparametric 
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mapping is performed, the whole element in X-Y plane is mapped onto an 
element in £“T ) plane, where the relative position of the predicted 
point and the closest corner does not change. A set of roots is solved 
for from the nonlinear system using this corresponding "closest corner" 
in £-t) plane as an Initial guess. The roots (£,i») are then found again 
using (0,0) the initial guess. If a root can be uniquely obtained 
inside the element by the initial guesses of (0,0) and the "closest 
corner”, then the point to be predicted is said to be located In a 
sufficiently undistorted part of the element. If the root obtained by 
initial guess of (0,0) is different from that by the closest corner, 
then it Implies there is an unacceptable geometric irregularity at the 
associated corner with reference to the physical coordinate plane. 
Considering element 4, the uniqueness of the roots seems valid for the 
guesses inside the element, but is ruined when the closest corner 
(-1,1) is used as a guess. This is because the corner (2. 5,2. 5) in 
physical coordinates, which is associated with the corner (-1,1) in 
natural coordinate, is collinear with the other 2 corners of the 
element. Similarly, the element 5 is also irregular at the corner 
(2. 5, 2. 5), so that the associated closest corner (-1,1) will lead to a 
different set of roots from that by the guess of (0,0). 

In the same example, if we try to predict the damage at (3. 2, 3. 2) 
just by the element 1, element 4 and element 5, we find not only the 
regular element 1 but also the irregular elements 4 and 5 can lead to a 
close approximation with the closest corners used as initial guesses, 
Table 3.4. This is because (3. 2, 3. 2) is much closer to these regular 
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element 

• 

roots 

• • 

roots 

• 

damage 

• • 

damage 

2345 

0.86,0.32 

0.86,0.32 

3. 87E+01 

3.87E+01 

1235 

0.76,0.36 

0.76,0.36 

3.84E+01 

3.84E+01 

1234 

0.78,0.31 

0.78,0.31 

3.84E+01 

3.84E+01 


note: 

root* : root (£,7)) obtained by initial guess of (0,0) 

root**: root (£,?)) obtained by initial guess of closest corner 
* * 
damage ; damage based on root 

•• *• 

damage : damage based on root 

exact damage : 39. 168 


Table 3.4 Predictions of the damage at (3. 2,3. 2) made by regular 
and irregular elements. 
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corners than to the distorted corners. 


3.2,2 Example 2 

This second example is presented to explain how to make a 
prediction when the predicted node is outside all possible elements. 
Suppose we are to predict the damage at (4,3) instead of (3,3), and the 
rest nodal data are the same as those in example 1. After finding the 
"scaled distances" and sorting, five different elements (see Figure 3.8 
a-e) are formed by taking any four distinct nodes from the Table 3-5. 


no. 

V 

t 

D 

scaled distance 
from ( v , tT— (4, 3) 

i 

4 

2 

38 

0.435 

2 

3 

4 

43 

0.591 

3 

2.5 

2.5 

20.625 

0.638 

4 

2 

1 

7 

1. 182 

5 

1 

2 

4 

1.276 


Table 3.5 The data in example 2 are sorted by scaled distances. 


Then by using the observation made on the results of example 1 that the 
predicted node should be inside the element if an initial guess of 
(0,0) is used for solving the nonlinear system such that the root is 
between -1 and 1, we find the node (4,3) is actually outside these 5 
elements (see Table 3.6 a-e). If linear interpolation functions are 
assumed applicable when the point to be predicted is outside of the 
element but is still sufficiently close to the element, then we can try 
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element 2 [1345] 


initial guess 

root 

damage 

0 , 0 

3.0 , 1.0 

3.72E+01 

1 . 1 

3.0 , 1.0 

3.72E+01 

1 ,-l * 

3.0 , 1.0 

3.72E+01 

-1 ,"1 

3.0 , 1.0 

3.72E+01 

-1 , 1 

3.0 , 1.0 

3.72E+01 



(b) 



element 3 [1245] 


initial guess 

root 

damage 

0 , 0 

1.37 , -0.26 

4.62E+01 

1 , 1 

1.37 , -0.26 

4. 62E+01 

1 ,-l * 

1.37 , -0.26 

4. 62E+01 

-1 ,-l 

1.37 , -0.26 

4. 62E+01 

-1 , 1 

1.37 , -0.26 

4.62E+01 



(c) 
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element 4 [1235] 

initial guess 

root 

damage 

■n 

1.44 , -0.14 

4.66E+01 


1.44 , -0.14 

4.66E+01 


1.44 , -0.14 

4.66E+01 


-1 , 3 

— 

watm 

-1 , 3 

— 


(d) 


element 5 [1234] 

initial guess 

root 

damage 

0 , 0 


4.66E+01 

1 , 1 

1.49 , -0.29 

4.66E+01 

1 ,-l * 

1.49 , -0.29 

4.66E+01 

-1 ,-l 

1.49 , -0.29 

4.66E+01 

-1 , 1 

diverge 

— 


(e) 

Table 3-6 Comparisons of the predicted damages in example 2, based on 
different elements and initial guesses (exact damage is 55) 
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to check the irregularities of the elements and the uniqueness of the 
mapping just as we treated the cases in example 1. We take the 4 
corners of each element in the natural coordinate plane as Initial 
guesses to solve the nonlinear systems, and then determine the 
associated damages (see Table 3.6 a-e). From these 5 tables, we find 
that both element 2 and element 3 with regular shapes, where the roots 
are uniquely determined, cam be used to obtain a prediction of the 
exact damage of 55. An approximation of D cam also be obtained from 
element 4 and element 5, though there are one or two corners that ruin 
the uniqueness of the the roots due to their geometric distortion. 

Like example 1, we find that the node (4,3) is closer to the regular 
corners of the elements 2 and 3 such that the uniqueness of the roots 
is still valid around these corners. Conversely, node (4,3) is closer 
to the irregular corner of the element 1, where 3 nodes are col linear 
and the root camnot converge. The corner that is closest to (4,3) with 
reference to natural coordinate plane is marked with ’ in each of the 
Table 3.6 a-e. A similar observation to that which was made 
considering example 1 can be made here. If the root obtained by the 
initial guess of the corner which is closest to the predicted node in 
the natural coordinate plame is the saune am that by the guess of (0,0), 
then the uniqueness of mapping is assumed satisfied and the 
approximation can be thought reasonable. Thus, except element 1, we 
have 4 possible approximations of D in this example. Because the 
predicton node should be as near the element we use as possible for 
applying linear interpolation functions, the prediction associated with 
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the root closest to the origin of the natural coordinate plane will be 
considered the most accurate. As shown in table 3.7, the distance 
associated with the 3rd element is the shortest and is marked with 
so the final approximation is determined on the basis of this element. 

3.3 Building a Criterion for Prediction 

Based on the 2 examples considered in this section, it is not 
difficult to extend the same approach to problems of higher DOF. 
Therefore, we can make a criterion for prediction by generalizing the 


element 

distance to 0 

2345 

— 

1345 

3. 162 

1245 

1.396 * 

1235 

1.451 

1234 

1.514 


Table 3.7 Comparison of the distances from each natural coordinate 
solved from the nonlinear solver to the origin 


observations made as follows : 

An ideal Isoparametric approximation model should be based on a least 
distorted element which contains the point to be predicted. However, 
it is not always possible to use undistorted elements in practical 
applications, especially for the cases of higher DOF. Accordingly, a 
unique isoparametric mapping over the whole element often cannot be 
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obtained. For this reason, only the uniqueness of the mapping over the 
relatively undistorted part of an element is used by comparing the 
roots from the nonlinear system when the initial guess of (0,0) is made 
and when the initial guess of the corner closest to the predicted node 
with reference to the natural coordinate are used. Thus, even for 
distorted elements, the "partial uniqueness" of the mapping can result 
in a reasonable function prediction as long as the predicted node 
either inside or outside the element is near an undistorted part of the 
element. When the node to be predicted is outside all possible 
elements, the final approximation will be determined by the root 
closest to the origin of the natural coordinate, if several possible 
approximations are available. 


104 



CHAPTER FOUR 


NUMERICAL TESTING AND RESULTS 


4.1 Sources of Data and the Testing Procedure 

This chapter will be devoted to testing the interpolation — 
extrapolation model using two sets of actual experimental data. The 
First set of data that is given Table 4. 1 was collected from 
experiments on debris Impact on the simulated bumper of the Space 
Station performed at Marshall Space Flight Center [22]. This data will 
be used for testing 3-DOF and 4-DOF interpolation models in the 
following sections. As shown in the Table 4.1, the independent 
parameters of t,d,0, and v represent the thickness of the bumper, the 
diameter of the debris particle, the Impact angle, and impact velocity, 
respectively. These parameters were illustrated in Figure 2.3. Based 
on these parameters, the dependent variables of the major diameter of 
the bumper hole (Dmj) and the minor diameter of the bumper hole (Dain) 
will be predicted. 

In addition, for evaluating the flexibility and the versatility 
of our model, a second set of data from the field of geology [23] will 
be used for testing the 5-DOF case, Table 4.2. The data in Table 4.2 
are related to the problem of determining the basin magnitude (Y), 
which essentially is a count of the number of sources in the basin, by 
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Table 4.1 : Listing of the experimental data of debris 


impact on the simulated bumper of space station 


t 

d 

0 

V 

Dmaj 

Drain 

0.063 

0.313 

45 

6.39 

0.81 

0.62 

0.04 

0.313 

45 

6.51 

0.76 

0.54 

0.04 

0.25 

45 

5.51 

0.53 

0.43 

0.04 

0.25 

45 

7.21 

0.53 

0.25 

0.04 

0.3 

65 

6.45 

0.94 

0.53 

0.04 

0.3 

65 

3.67 

0.87 

0.47 

0.04 

0.35 

65 

3.04 

0.98 

0 . 49 

0.063 

0.25 

45 

4 . 11 

0.61 

0.48 

0.063 

0.25 

45 

4.59 

0.65 

0.49 

0.063 

0.25 

45 

5.30 

0.6 

0.5 

0.063 

0.3 

65 

5.85 

0.88 

0.59 

0.063 

0.3 

65 

7.08 

1.02 

0.64 

0.063 

0.25 

65 

6.4 

0.87 

0.52 

0.063 

0.25 

65 

7.4 

0.77 

0.57 

0.063 

0.35 

65 

5.7 

1 . 13 

0.66 

0.063 

0.35 

65 

6.8 

1.4 

0.68 

0.063 

0.35 

45 

5.85 

0.86 

0.68 

0.04 

0.35 

45 

6.4 

0.8 

0.6 

0.04 

0.35 

45 

6.76 

0.84 

0.59 

0.04 

0 . 187 

45 

6.36 

0.48 

0 . 4 

0.04 

0 . 187 

45 

5.89 

0.54 

0.42 

0.04 

0 . 187 

45 

4.57 

0.45 

0.36 

0.032 

0.25 

45 

5.52 

0.6 

0.46 

0.032 

0.25 

45 

7 . 12 

0.6 

0 . 43 

0.063 

0 . 187 

45 

4.41 

0.48 

0.39 

0.063 

0 . 187 

45 

3.23 

0.47 

0.36 

0.063 

0.25 

45 

2.95 

0.54 

0.43 

0.063 

0.313 

45 

4.59 

0.73 

0 . 58 

0.063 

0.313 

45 

4.34 

0.71 

0.57 

0.04 

0 . 187 

45 

2.88 

0.4 

0.3 

0 . 125 

0 . 187 

45 

3.61 

0.4 

0.31 

0.04 

0.25 

45 

4.47 

0.53 

0.42 

0.04 

0.313 

45 

7.0 

0.77 

0.52 

0.04 

0.313 

45 

6.98 

0.84 

0.57 

0.04 

0.375 

45 

6.49 

0.83 

0.69 
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Table 4.2 : Listing of geological observed data 


used to determine basin magnitude 


elevation 

rel ief 

720 

570 

670 

610 

860 

550 

870 

610 

730 

570 

690 

590 

880 

640 

760 

690 

820 

600 

720 

480 

670 

670 

660 

600 

830 

660 

780 

620 

750 

740 

770 

630 

750 

570 

750 

580 

740 

760 

750 

740 

750 

760 

740 

770 

940 

510 

700 

600 

810 

580 

920 

500 

920 

490 

790 

605 

860 

550 

860 

630 

880 

520 

780 

460 

720 

440 

780 

300 

700 

460 

680 

520 

820 

520 


length of 
area the stream 


7 

154 

3 

80 

11 

84 

11 

122 

14 

185 

12 

200 

11 

170 

28 

340 

5 

100 

3 

80 

19 

290 

5 

90 

18 

260 

17 

111 

15 

184 

21 

227 

4 

60 

20 

259 

9 

62 

6 

95 

11 

105 

32 

350 

21 

232 

23 

266 

44 

390 

13 

142 

12 

145 

33 

253 

23 

241 

87 

702 

37 

288 

17 

162 

8 

67 

3 

52 

10 

121 

26 

220 

8 

123 


(continued) 


drainage basin 
density magnitude 


2200 

14 

2667 

6 

763 

5 

1110 

7 

1321 

11 

1667 

14 

1545 

12 

1215 

18 

2000 

6 

2667 

5 

1526 

17 

1800 

5 

1444 

22 

652 

7 

1227 

15 

1080 

17 

1500 

5 

1295 

18 

689 

14 

1583 

21 

954 

22 

1094 

23 

1105 

28 

1156 

42 

886 

22 

1092 

10 

1208 

11 

766 

12 

1048 

13 

807 

31 

778 

18 

953 

13 

838 

4 

1733 

5 

1210 

9 

846 

13 

1537 

10 


107 


Table 4.2 (continued) 


length of drainage basin 


levat Ion 

relief 

area 

the stream 

density 

magnitude 

710 

520 

24 

238 

992 

13 

800 

440 

19 

231 

1216 

13 

700 

510 

16 

178 

1113 

11 

675 

570 

18 

168 

933 

12 

740 

510 

8 

65 

812 

4 

740 

520 

31 

334 

1078 

17 

770 

600 

21 

184 

876 

9 

820 

520 

11 

136 

1237 

8 

850 

490 

22 

233 

1059 

13 

820 

629 

34 

410 

1206 

22 

820 

510 

11 

149 

1354 

10 

680 

640 

46 

348 

757 

19 

660 

789 

55 

382 

695 

27 
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considering the 6 following Independent variables; the elevation of the 
basin outlet (x ), the relief of the basin (x ), the basin area (x ), 

1 2 3 

total length of the stream in the basin (x^), and drainage density 

(x ), which is defined as total length of the streams in the basin 
5 

divided by the basin area. 

The interpolation/extrapolation procedure was tested in following 
manner. The dependent variable (say, diameter or basin magnitude) of 
each set of data is assumed unknown and then is predicted based on 
the remaining data. For example, in Table 4-1, suppose the major 
diameter (or minor diameter) associated with the first set of data is 
set as unknown, then the rest of 34 sets of data will be used to predict 
it. Similarly, the major diameters associated with the second, ..., etc. 
set of the data will be predicted by the others. 

For simplicity, each prediction was based on the 2 n + 1 closest 
elements ( n = number of DOF ). As was suggested in chapter 3, the 
final prediction was based on the element that could produce unique 
natural coordinate roots after solving the nonlinear mapping equations 
by two initial guesses. If more than one element produced unique roots 
in the test, then the final prediction was taken as that produced by 
the element whose natural coordinate roots were closest to the origin 
of the natural coordinate system. For cases where no unique mapping 
was found in all the tested elements, then the roots closest to the 
origin were used to make the prediction. 

Predicted values will now be compared with that measured. A 
computer program ’ ISOMODEL ’ that implements the algorithms derived in 
Chapters 2 and 3 was written in Turbo Pascal 5.5 to perform all the 
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tests on IBM PC 386 machine. This program Is listed in the appendix. 
Because multiple linear regression is a standard and well accepted 
technique for this sort of problem, the predicted results generated by 
the statistical software package SAS [24] using the IBM 3090 are 
compared with the results generated by the Isoparametric model. The 
linear functions that SAS used are of the form : 

Y = a + aX + aX + ... + aX 

0 11 2 2 n n 

where Y is dependent variable, X^ (i * l..n) are Independent variables, 
and a^ (i = 0..n) are coefficients to be determined by the Least Square 
method. 

4.2 Testing for 3-D Case and the Results 

In this section, the 3-D0F case is tested by using a subset of 
Table 4.1, where there are just 26 sets of data that correspond to 0 = 
45. We predict the major diameter (Dmaj) first, and then the minor 
diameter (Dmin) by the same parameters. The results of prediction on 
major diameter and minor diameter, compared with those predicted by 
SAS, are respectively shown in Table 4.3 and 4.4. 

As Table 4.3 and Table 4.4 show, the average error of the 
isoparametric model (9.37%) is greater than that of SAS (5.61%) In 
absolute value for the prediction of the major diameter. For the 
prediction of the minor diameter, SAS also has less average error 
(9.64%) than the isoparametric model (11.68%). The few wild predicted 
values of the isoparametric model can be attributed to a fatal 
distortion of elements or possibly to some scatter in the experimental 
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Table 4.3 : Listing of the prediction of major diameter 
by the isoparametric model and SAS vs. the 
measured value 


No. 

MEASURED 
Dma J 

PREDICTED 

IS0M0DEL 

BY 

ERROR?. 

PREDICTED BY 
SAS 

ERR0R% 

i 

0.81 

0.79 

-2.47 

0.76 

-6. 17 

2 

0.76 

0.76 

0.00 

0.75 

-1.32 

3 

0.53 

0.58 

11.32 

0.60 

13.21 

4 

0.53 

0.67 

26.42 

0.63 

18.87 

5 

0.61 

0.62 

1.64 

0.58 

-4.92 

6 

0.65 

0.62 

-4.62 

0.59 

-9.23 

7 

0.6 

0.68 

13.33 

0.60 

0.00 

8 

0.86 

0.86 

0.00 

0.83 

-3.49 

9 

0.8 

0.82 

2.50 

0.83 

3.75 

10 

0.84 

0.86 

2.38 

0.83 

-1.19 

11 

0.48 

0.57 

18.75 

0.48 

0.00 

12 

0.54 

0.47 

-12.96 

0.47 

-12.96 

13 

0.45 

0.48 

-2.22 

0.45 

0.00 

14 

0.6 

0.53 

-11.67 

0.60 

0.00 

15 

0.6 

0.48 

-20.00 

0.62 

3. 33 

16 

0.48 

0.56 

16.67 

0.45 

-6.25 

17 

0.47 

0.46 

-2. 13 

0.43 

-8. 51 

18 

0.54 

0.51 

-5.56 

0.56 

3.70 

19 

0.73 

0.72 

0.00 

0.72 

-1.37 

20 

0.71 

0.72 

1.41 

0.72 

1.41 

21 

0.4 

0.37 

-7.50 

0.42 

5.00 

22 

0.4 

0.47 

17.50 

0.46 

15.00 

23 

0.53 

0.29 

-45.28 

0.58 

9.43 

24 

0.77 

0.84 

9.09 

0.76 

-1.30 

25 

0.84 

0.77 

-8.33 

0.76 

-9.52 

26 

0.83 

0.88 

0.00 

0.88 

6.02 




Ave= 9.37% 


Ave=5. 61% 
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Table 4.4 : Listing of the prediction of minor diameter 
by the isoparametric model and SAS vs. the 
measured value 


No. 

MEASURED 
Dma J 

PREDICTED 

IS0M0DEL 

ERR0R% 

PREDICTED BY 
SAS 

ERROR 0 /. 

1 

0.62 

0.64 

3.23 

0.56 

-9.68 

2 

0.54 

0.56 

3.70 

0.55 

1.85 

3 

0.43 

0.45 

4.65 

0.45 

4.65 

4 

0.25 

0.51 

104.00 

0.44 

76.00 

5 

0.48 

0.47 

-2.08 

0.46 

-4. 17 

6 

0.49 

0.49 

0.00 

0.46 

-6. 12 

7 

0.50 

0.51 

2.00 

0.46 

-8.00 

8 

0.68 

0.66 

-2.94 

0.63 

-7.35 

9 

0.60 

0.65 

8.33 

0.62 

3.33 

10 

0. 59 

0.62 

5.08 

0.62 

5.08 

11 

0.40 

0.44 

10.00 

0.33 

-17.50 

12 

0.42 

0.39 

-7. 14 

0.34 

-19.05 

13 

0.36 

0.37 

2.78 

0.34 

-5. 56 

14 

0. 46 

0. 43 

-6.52 

0.44 

-4.35 

15 

0.43 

0. 16 

-62.79 

0.43 

0.00 

16 

0.39 

0.40 

2.56 

0.35 

-10.26 

17 

0.36 

0.35 

-2.78 

0.36 

0.00 

18 

0.43 

0.46 

6.98 

0.47 

9.30 

19 

0.58 

0.58 

0.00 

0.57 

-1.72 

20 

0.57 

0.57 

0.00 

0.57 

0.00 

21 

0.30 

0.30 

0.00 

0.35 

16.67 

22 

0. 31 

0.36 

16. 13 

0.37 

19.35 

23 

0.42 

0.33 

-21.43 

0.45 

7. 14 

24 

0.52 

0.57 

9.62 

0.55 

5.77 

25 

0.57 

0.52 

-8.77 

0.55 

-3.51 

26 

0.69 

0.62 

-10. 14 

0.66 

-4.35 


Ave= 11 . 68% Ave=9 . 64% 
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data. The comparisons of the results are shown in Figure 4. 1 and 
Figure 4.2. 

4.3 Testing for 4-D Case and the Results 

For the 4-DOF case, we take all the 4 independent parameters and 
all the 35 sets of data in Table 4.1 into consideration to test the 
isoparametric model. The results are listed in Table 4.5 and Table 
4.6. By comparing the average errors, the isoparametric model has less 
accurate predictions them SAS in major diameter but more accurate than 
SAS in minor diameter. Figure 4.3 and 4.4 also shows that the 
predictions made by SAS are more uniformly scattered along the 45 
degree line than those by the isoparametric model. The few wild 
predicted values of the isoparametric model could be caused by serious 
distortion of elements as well as the numerical errors from the 
nonlinear system solver due to the higher order nonlinear terms when 
Newton’s method applied. The inconsistency of part of the experimental 
data could also exaggerate the deviations for both models. 
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Figure 4.1 Predictions by the isoparametric model and SAS versu 
observed major diameter for the 3-DOF test. Diagonal 
line indicates predicted equal observed. 
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line indicates predicted 



Table 4.5 : Listing of the prediction of major diameter 
by the Isoparametric model and SAS vs. the 
measured value 


No. 

MEASURED 
Dma J 

PREDICTED 

IS0M0DEL 

ERR0R7. 

PREDICTED 

SAS 

BY 

ERR0R% 

1 

0.81 

0.80 

-1.23 

0.78 

-3. 70 

2 

0.76 

0.73 

-3.95 

0.76 

0.00 

3 

0.53 

0.59 

11.32 

0.59 

11.32 

4 

0.53 

0.75 

41.51 

0.63 

18.87 

5 

0.94 

0.76 

-19. 15 

0.97 

3. 19 

6 

0.87 

0.90 

3.45 

0.91 

4. 60 

7 

0.98 

0.95 

-3.06 

1.01 

3.06 

8 

0.61 

0.62 

1.64 

0.58 

-4.92 

9 

0.65 

0.61 

-6. 15 

0.59 

-9.23 

10 

0. 60 

0.70 

16.67 

0.61 

1.67 

11 

0.88 

0. 98 

11.36 

0.98 

11. 36 

12 

1.02 

1.08 

5.88 

1.01 

-0.98 

13 

0.87 

0.74 

-14.94 

0.87 

0.00 

14 

0.77 

1. 17 

51.95 

0.89 

15. 58 

15 

1. 13 

0.94 

-16.81 

1.09 

-3.54 

16 

1.40 

1.05 

-25. 00 

1. 12 

-20.00 

17 

0.86 

0.99 

15. 12 

0.85 

-1. 16 

18 

0.80 

0.73 

-8.75 

0.84 

5.00 

19 

0.84 

0.82 

-2.38 

0.85 

1. 19 

20 

0. 48 

0.54 

12.50 

0.46 

-4. 17 

21 

0.54 

0.49 

-9.26 

0.45 

-16.67 

22 

0.45 

0.48 

6.67 

0.42 

-6.67 

23 

0.60 

0.58 

-3.33 

0.58 

-3.33 

24 

0.60 

0.76 

26.67 

0.62 

3.33 

25 

0.48 

0.63 

31.25 

0.44 

-8.33 

26 

0.47 

0.48 

2. 13 

0.41 

-12.77 

27 

0.54 

0.61 

12.96 

0.55 

1.85 

28 

0.73 

0.73 

0.00 

0.74 

1.37 

29 

0.71 

0.71 

0.00 

0.73 

2.82 

30 

0.4 

0.45 

12.50 

0.38 

-5.00 

31 

0.4 

0.77 

92.50 

0.48 

20.00 

32 

0.53 

0.61 

15.09 

0.57 

7.55 

33 

0.77 

0.84 

9.09 

0.77 

0.00 

34 

0.84 

0.77 

-8.33 

0.77 

-8.33 

35 

0.83 

0.83 

0.00 

0.90 

8.43 



Ave=14. 367. 


Ave=6. 57% 
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Table 4.6 : Listing of the prediction of minor diameter 
by the Isoparametric model and SAS vs. the 
measured value 


No. 

MEASURED 
Dma j 

PREDICTED 

IS0M0DEL 

ERROR/. 

PREDICTED 

SAS 

BY 

ERROR/. 

i 

0.62 

0.66 

6.45 

0.58 

-6. 45 

2 

0.54 

0.55 

1.85 

0.56 

3.70 

3 

0. 43 

0.45 

4.65 

0.44 

2.33 

4 

0.25 

0.48 

92.00 

0.47 

88.00 

5 

0.53 

0.54 

1.89 

0.56 

5.66 

6 

0.47 

0.46 

-2. 13 

0.51 

8.51 

7 

0.49 

0.52 

6. 12 

0.57 

16.33 

8 

0.48 

0.47 

-2.08 

0.45 

-6.25 

9 

0. 49 

0.49 

0.00 

0.46 

-6. 12 

10 

0.50 

0.52 

4.00 

0.47 

-6.00 

11 

0.59 

0.59 

0.00 

0.57 

-3.39 

12 

0.64 

0.62 

-3.13 

0.60 

-6.25 

13 

0.52 

0.55 

5.77 

0.51 

-1.92 

14 

0.57 

0.54 

-5.26 

0.53 

-7.02 

15 

0.66 

0.65 

-1.52 

0.65 

-1.52 

16 

0.68 

0.70 

2.94 

0.67 

-1.47 

17 

0.68 

0.66 

-2.94 

0.63 

-7.35 

18 

0.60 

0.61 

1.67 

0.61 

1.67 

19 

0.59 

0.60 

1.70 

0.62 

5.08 

20 

0.40 

0.42 

5.00 

0.36 

-10.00 

21 

0.42 

0.41 

-2.38 

0.36 

-14.29 

22 

0.36 

0.37 

2.78 

0.33 

-8.33 

23 

0.46 

0.46 

0.00 

0.43 

-6.52 

24 

0.43 

0.47 

9.30 

0.46 

6.98 

25 

0.39 

0.40 

2.56 

0.36 

-7.69 

26 

0.36 

0.39 

8.33 

0.34 

-5.56 

27 

0.43 

0.48 

11.63 

0.43 

0.00 

28 

0.58 

0.58 

0.00 

0.55 

-5. 17 

29 

0.57 

0.57 

0.00 

0.55 

-3.51 

30 

0.3 

0.36 

20.00 

0.30 

0.00 

31 

0.31 

0.46 

48.39 

0.43 

38.71 

32 

0.42 

0.40 

-4.76 

0.43 

2.38 

33 

0.52 

0.57 

9.62 

0.56 

7.69 

34 

0.57 

0.52 

-8.77 

0.56 

-1.75 

35 

0.69 

0.64 

-7.25 

0.65 

-5.80 




Ave=8. 20% 


Ave=8.84*/. 
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Figure 4.3 Predictions by the isoparametric model and SAS versu 
observed major diameter for the 4-DOF test. Diagonal 
line indicates predicted equal observed. 
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4.4 Testing for 5-D Case 

In this section, the 50 sets of 5-DOF data of Table 4.2 were used 
to test the isoparametric model. The results are listed in Table 4.7. 
The isoparametric model and SAS have much higher average errors (48.20% 
and 25.48%, respectively) than in the previous cases. Referring to 
Figure 4.5, it is can be seen that SAS can make a closer prediction to 
the observed values than the isoparametric model. More wild predicted 
values by the isoparametric model appear than those in 3-DOF and 4-DOF 
cases. This could be partially due to the effects of distortion on 
elements with higher DOF, and partially due to the experimental 
scatter. 
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Table 4.7 : Listing of the prediction of basin magnitude 
by the isoparametric model and SAS vs. the 
measured value 


MEASURED PREDICTED PREDICTED BY 


No. 

Dma J 

ISOMODEL 

ERROR?. 

SAS 

ERROR?. 

1 

14 

10. 17 

- 27.36 

11.53 

-17.64 

2 

6 

4.04 

-32. 67 

7.43 

23.83 

3 

5 

5.27 

5.40 

8.69 

73.80 

4 

7 

-2. 15 

-130.71 

11.76 

68.00 

5 

11 

11.53 

4.82 

13.81 

25.55 

6 

14 

11. 19 

-20.07 

14.59 

6.79 

7 

12 

15.81 

31.75 

14.65 

22.08 

8 

18 

17.30 

-3.89 

23.57 

30.94 

9 

6 

4.22 

-29.67 

9.57 

59.50 

10 

5 

2.03 

-59.40 

4.63 

-7.40 

11 

17 

11.27 

-33.71 

21.26 

25.06 

12 

5 

6.52 

25.00 

9.26 

85.20 

13 

22 

17.34 

-21. 18 

19.56 

-11.09 

14 

7 

-6.4 

-191.43 

10.76 

53.71 

15 

15 

36.67 

144. 47 

17.42 

16. 13 

16 

17 

23.50 

38.24 

16.76 

-1.41 

17 

5 

7.92 

58.40 

7.56 

51.20 

18 

18 

26.86 

49.22 

17.56 

-2.44 

19 

14 

20.58 

47.00 

12.37 

-11.64 

20 

21 

12.60 

-40.00 

12.91 

-38.52 

21 

22 

14.29 

-35.05 

14. 16 

-35.64 

22 

23 

22.54 

-2.00 

25.27 

9.87 

23 

28 

12.78 

-54. 36 

14.52 

-48. 14 

24 

42 

15.94 

-62.05 

18.01 

-57. 12 

25 

22 

19.25 

-12.50 

21.41 

-2.68 

26 

10 

10.02 

0.20 

10.28 

2.80 

27 

11 

8.29 

-24. 64 

10.25 

-6.81 

28 

12 

0.90 

-92.50 

15.77 

31.42 

29 

13 

13.53 

4.08 

15.59 

19.92 

30 

31 

51.63 

66.55 

33.06 

6.65 

31 

18 

5.34 

-70.33 

15.30 

-15.00 

32 

13 

10.73 

-17.46 

10.00 

-23.08 

33 

4 

10.44 

161.00 

5.65 

41.25 

34 

5 

-1.3 

-126.00 

0.91 

-81.80 

35 

9 

4.0 

-55. 56 

8.37 

-7.00 


(continued) 
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(Table 4.7 continued) 


MEASURED PREDICTED PREDICTED BY 


No. 

Dma J 

ISOMODEL 

ERROR?. 

SAS 

ERROR’/. 

36 

13 

-3.8 

-129.23 

13. 18 

1. 38 

37 

10 

7.02 

-29.80 

9.64 

-3.60 

38 

13 

14.01 

7.77 

14.55 

11.92 

39 

13 

18.77 

44.38 

13. 10 

0.77 

40 

11 

9.56 

-13.09 

12.00 

9.09 

41 

12 

12.43 

3.58 

12.54 

4. 50 

42 

4 

6.23 

55.75 

7. 10 

77.50 

43 

17 

23.93 

40.76 

19.05 

12.06 

44 

9 

10.04 

11.56 

13.71 

52.33 

45 

8 

10.72 

34.00 

10.43 

30.38 

46 

13 

29.29 

125.31 

13.97 

7.46 

47 

22 

0. 63 

-97. 13 

25.50 

15. 91 

48 

10 

8.80 

-12.00 

10.82 

8.20 

49 

19 

21. 12 

11. 16 

19.67 

3. 53 

50 

27 

22.76 

-15.70 

23. 13 

-14.33 




Ave=48.20% 


Ave=25 . 48% 
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Figure 4.5 Predictions by the isoparametric model and SAS versus 
observed basin number for the 5-DOF test. Diagonal 
line indicates predicted equal observed. 



CHAPTER FIVE 


CONCLUSIONS AND RECOMMENDATIONS 


This thesis has proposed a unique technique for multivariable 
analysis based on extending Isoparametric concept of FBI. By extending 
the concept of shape functions defined by geometric coordinates In 2-D 
and 3-D space to that of those defined by physical coordinates In n-D 
space, a multivariable Interpolation model was developed to analyze 
scattered data with arbitrary DOF. A Pascal program called ISOMODEL, 
based on the methodology derived In Chapter 3, was tested for the 3-D, 
4-D and 5-D cases by two sets of actual experimental data In Chapter 4. 
Compared with the results from the statistical software SAS, the 
isoparametric model was less accurate In the 3-D and 4-D cases, but the 
difference of the average errors of these two models is not much. 
However, In the cases of higher DOF, the Isoparametric model Is far 
less accurate than SAS. Geometrically, this is because the 
isoparametric model Is too sensitive to the distortion of the elements. 
Technically, the Isoparametric model Involves solving nonlinear systems 
of equations, which often causes remarkable numerical errors in 
treating higher order nonlinear terms. At the same time, errors 
originating from the scatter in experimental observations and 
measurements could also Increase the prediction errors. 
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It should be noticed that the points to be predicted were 
excluded and then estimated by the remaining sets of data when the 
Isoparametric model was tested. On the contrary, the points to be 
predicted, together with the remaining data points, were used to 
determine the linear regression model when SAS was tested. This could 
be another factor to cause some bias in predictions between the 
Isoparametric model and SAS. 

In order to improve the weak points of the isoparametric model, 
further study is required. It is recommended that an alternate way be 
developed for finding the natural coordinates in the isoparametric 
mapping such that the errors of solving nonlinear systems can be 
reduced. Finding a more effective way to find the least distorted 
elements from the nodal points could help lower the possibility of the 
wild predictions. 
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Appendix 2 A Listing of the Computer Code. (Listing not provided if floppy 

DISKS INCLUDED) 



mxmiwim mm 
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PRECEDING PAGE BLANK NOT FILMED 
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MLIBLAST.BAS 


MLIBLAST . BAS 

Source code for the main program that runs the other programs 


CLS 

COLOR 15, 0 
LOCATE 12, 37 
PRINT "MLIBLAST" 

seed% = ((TIMER * 65536) / 86400) - 32768 
RANDOMIZE seed% 

FOR 1% = 1 TO 2000 
testcolor# = RND 

IF testcolor# >= .666 THEN colornum% = 15 
IF testcolor# <= .333 THEN colomum% = 9 

IF testcolor# > .333 AND testcolor# < .666 THEN colornum% = 12 
COLOR colornum% , 0 
row% = 1! + RND * 221 
col% = 11 + RND * 791 

IF row% < 11 OR row% >13 THEN 
LOCATE row%, col% 

PRINT CHR$ ( 219 ) 

ELSE 

IF col% < 36 OR col% > 45 THEN 
LOCATE row%, col% 

PRINT CHR$ ( 219 ) 

END IF 
END IF 
NEXT 1% 

MainMenu : 

COLOR 15, 9 
CLS 

LOCATE 1, 3 

PRINT "Please Enter The "; 

COLOR 12, 9 
PRINT "Number "? 

COLOR 15, 9 

PRINT "Associated With the Desired Actions" 

LOCATE 5, 3 
COLOR 12, 9 
PRINT " 1 . " ; 

COLOR 15, 9 

PRINT "Add Data To The Database" 

LOCATE 7, 3 
COLOR 12, 9 
PRINT " 2 . " ; 

COLOR 15, 9 

PRINT "Predict MLI Damage" 

LOCATE 9, 3 
COLOR 12, 9 
PRINT " 3 . " ? 

COLOR 15, 9 
PRINT "Quit" 

LOCATE 13, 3 
COLOR 12, 9 
INPUT Choice% 

IF Choice% < 1 OR Choice% > 3 THEN GOTO MainMenu 
/ 

IF Choice% = 1 THEN 


J PBHMWB mm 


PRSC£PLNG PAGE BLANK NOT F"LMP*L» 



MLIBLAST.BAS 

SHELL " database.exe" 

GOTO MainMenu 
END IF 

f 

IF Choice% » 2 THEN 
CLS 

COLOR 15, 9 
LOCATE 1, 3 

PRINT "Please Enter The " ; 

COLOR 12, 9 
PRINT "Number ' ; 

COLOR 15, 9 

PRINT "Associated With the Desired Prediction Scheme:" 
LOCATE 5, 3 
COLOR 12, 9 
PRINT " 1 . " ; 

COLOR 15, 9 

PRINT "1/R~N Interpolation Scheme" 

LOCATE 7, 3 
COLOR 12, 9 
PRINT " 2 . " ; 

COLOR 15, 9 

PRINT "Polynomial Interpolation Scheme" 

LOCATE 9, 3 
COLOR 12, 9 
PRINT " 3 . " ; 

COLOR 15, 9 

PRINT "Nondimensional Parameter Scheme" 

LOCATE 11, 3 
COLOR 12, 9 
PRINT " 4 . " ; 

COLOR 15, 9 

PRINT "Return to the Main Menu" 

LOCATE 13, 3 
COLOR 12, 9 
INPUT Choice% 

IF Choice% = 1 THEN SHELL "invrmeth.exe" 

IF Choice% = 2 THEN SHELL "polymeth.exe" 

IF Choice% - 3 THEN SHELL "nondimen.exe" 

IF Choice% = 4 THEN GOTO MainMenu 

GOTO MainMenu 
END IF 

IF Choice% = 3 THEN 
COLOR 15, 0 
CLS 
END 
END IF 
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DATABASE.BAS 


DATABASE. BAS 

Source code for the data base program 


ON KEY (10) GOSUB DatalnputProblem 
AddData : 

COLOR 11, 0 
CLS 

LOCATE 1, 3 

PRINT "Enter MLI Test Data File Name: 
COLOR 12, 0 

INPUT " MLITestDataFile$ 

OPEN MLITestDataFile$ FOR APPEND AS #1 
AddData 1 : 

KEY (10) ON 

GOSUB DisplayTemplate 

9 

'prompt user for data 


COLOR 12, 

0 


LOCATE 25 

, 1 


PRINT "Please Enter Data At The Cursor 

- Press F10 and ENTER To Restart Input" 

COLOR 12, 

9 

0 


LOCATE 2, 

10 


INPUT , 

9 

TestID$ 


LOCATE 2, 

37 


INPUT "", 

9 

DataSource$ 


9 

LOCATE 2, 

62 


INPUT " ", 

9 

TestDateS 


9 

LOCATE 5, 

18 


INPUT "", 

9 

BumperMaterial$ 


9 

LOCATE 5, 

48 


INPUT " ", 

f 

BumperThickness# 


9 

LOCATE 5, 

74 


INPUT "", 

9 

BumperS tandOf f # 


9 

LOCATE 8, 

25 


INPUT "", 

9 

PressureWallMaterial$ 


9 

LOCATE 8, 

74 


INPUT "", 

9 

PressureWallThickness# 


9 

LOCATE 11 

, 22 


INPUT "", 

Pro jectileMaterial$ 
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DATABASE.BAS 


LOCATE 11, 74 

INPUT ProjectileDiameter# 


LOCATE 14, 15 

INPUT ImpactAngle# 


LOCATE 14, 70 

INPUT Pro jectileVelocity# 


LOCATE 17, 33 

INPUT BumperMa jorAxis# 


LOCATE 17, 70 

INPUT BumperMinorAxis# 


LOCATE 20, 20 

INPUT MLIHoleDiam# 


LOCATE 20, 70 

INPUT MLIMassLoss# 


LOCATE 23, 39 

INPUT PressWallMajAxis# 


LOCATE 23, 70 

INPUT PressWallMinAxis# 


KEY ( 10 ) OFF 
LOCATE 25, 1 
PRINT " 

LOCATE 25, 1 
COLOR 9, 0 

PRINT "OK to write this data to the database (y/n)7 
COLOR 12, 0 
INPUT " " , answer$ 

COLOR 15, 0 

answer$ * LCASE$ (answer$ ) 

IF answer $ = "y" THEN 
GOSUB WriteDataToFile 
CLS 

LOCATE 25, 1 
COLOR 9, 0 

PRINT "Do you wish to enter more data at this time (y/n)? 
COLOR 12, 0 
INPUT " " , answer$ 
answer$ = LCASE$(answer$) 

COLOR 15, 0 

IF answer$ = B y M THEN 132 

GOTO AddDatal 
ELSE 



DATABASE.BAS 


GOTO Finish 
END IF 
ELSE 
CLS 

LOCATE 25, 1 
COLOR 9, 0 

PRINT "Do you wish to enter more data at this time (y/n)? 
COLOR 12, 0 
INPUT " " , answer$ 

COLOR 15, 0 
IF answers = "y" THEN 
GOTO AddDatal 
ELSE 

GOTO Finish 
END IF 
END IF 
Finish: 

END 

DisplayTemplate : 

CLS 

COLOR 11, 0 

r 

'field 1 

r 

LOCATE 2, 1 
PRINT "Test ID" 

BoxRow% = 1 
BoxColumn% = 9 
BoxLength% = 10 
GOSUB BoxDraw 

r 

'field 2 
/ 

LOCATE 2, 24 
PRINT "Data Source" 

BoxRow% = 1 
BoxColumn% = 36 
BoxLength% = 10 
GOSUB BoxDraw 

r 

'field 3 
/ 

LOCATE 2, 51 
PRINT "Test Date" 

BoxRow% = 1 
BoxColumn% = 61 
BoxLength% = 18 
GOSUB BoxDraw 

t 

'field 4 
# 

LOCATE 5, 1 

PRINT "Bumper Material" 

LOCATE 6, 7: PRINT "Name" 

BoxRow% ■ 4 
BoxColumn% = 17 
BoxLength% = 10 

GOSUB BoxDraw 133 

f 


'field 5 



DATABASE. BAS 


LOCATE 5/ 30 

PRINT "Bumper Thickness" 

LOCATE 6, 36 
PRINT "(in)" 

BoxRow% = 4 
BoxColumn% =47 
BoxLength% = 6 
GOSUB BoxDraw 

f 

'field 6 

f 

LOCATE 5, 56 

PRINT "Bumper Stand-Off" 

LOCATE 6, 62 
PRINT "(in)" 

BoxRow% = 4 
BoxColumn% = 73 
BoxLength% = 6 
GOSUB BoxDraw 

r 

'field 7 
/ 

LOCATE 8, 1 

PRINT "Pressure Wall Material" 

LOCATE 9, 10s PRINT "Name" 

BoxRow% = 7 
BoxColumn% = 24 
BoxLength% = 10 
GOSUB BoxDraw 

f 

'field 8 
/ 

LOCATE 8, 48 

PRINT "Pressure Wall Thickness" 

LOCATE 9, 58 
PRINT " ( in ) " 

BoxRow% = 7 
BoxColumn% =73 
BoxLength% = 6 
GOSUB BoxDraw 

r 

' field 9 

r 

LOCATE 11, 1 

PRINT "Projectile Material" 

LOCATE 12, 9: PRINT "Name" 

BoxRow% = 10 
BoxColumn% = 21 
BoxLength% = 10 

GOSUB BoxDraw 
1 

'field 10 

f 

LOCATE 11, 53 

PRINT "Projectile Diameter" 

LOCATE 12, 61 
PRINT "(in)" 

BoxRow% =10 134 

BoxColumn% =73 



DATABASE.BAS 

BoxLength% = 6 
GOSUB BoxDraw 

9 

'field 10 

9 

LOCATE 14, 1 

PRINT "Impact Angle" 

LOCATE 15, 2 
PRINT " (degrees ) " 

BoxRow% = 13 
BoxColumn% = 14 
BoxLength% = 6 
GOSUB BoxDraw 

t 

'field 11 

9 

LOCATE 14, 49 

PRINT "Projectile Velocity" 
LOCATE 15, 55 
PRINT "(km/sec)" 

BoxRow% = 13 
BoxColumn% = 69 
BoxLength% = 10 
GOSUB BoxDraw 

9 

'field 12 

9 

LOCATE 17, 1 

PRINT "Bumper Hole Size -" 
LOCATE 17, 21 
PRINT "Major Axis" 

LOCATE 18, 24 
PRINT "(in)" 

BoxRow% = 16 
BoxColumn% = 32 
BoxLength% = 10 
GOSUB BoxDraw 

t 

'field 13 

f 

LOCATE 17, 58 
PRINT "Minor Axis" 

LOCATE 18, 61 
PRINT "(in)" 

BoxRow% = 16 
BoxColumn% = 69 
BoxLength% = 10 
GOSUB BoxDraw 

9 

'field 14 

9 

LOCATE 20, 1 

PRINT "MLI Hole Diameter" 
LOCATE 21, 8 
PRINT "(in)" 

BoxRow% = 19 
BoxColumn% = 19 
BoxLength% = 10 
GOSUB BoxDraw 



DATABASE.BAS 

'field 15 
/ 

LOCATE 20, 55 
PRINT "MLI Mass Loss" 

LOCATE 21, 58 
PRINT "(grams)" 

BoxRow% = 19 
BoxColumn% = 69 
BoxLength% = 10 
GOSUB BoxDraw 

9 

'field 16 

9 

LOCATE 23, 1 

PRINT "Pressure Wall Hole Size - Major Axis" 
LOCATE 24, 30 
PRINT "(in)"; 

BoxRow% = 22 
BoxColumn% = 38 
BoxLength% = 10 
GOSUB BoxDraw 

r 

'field 17 

f 

LOCATE 23, 58 
PRINT "Minor Axis" 

LOCATE 24, 61 
PRINT "(in)"; 

BoxRow% = 22 
BoxColumn% = 69 
BoxLength% =10 
GOSUB BoxDraw 
RETURN 
BoxDraw: 

COLOR 9, 0 

LOCATE BoxRow%, BoxColumn% 

BoxLine$ = CHR$(201) 

FOR i% = 1 TO BoxLength% 

BoxLine$ = BoxLine$ + CHR$(205) 

NEXT i% 

BoxLine$ ■ BoxLine$ + CHR$(187) 

PRINT BoxLine$; 

LOCATE BoxRow% + 1, BoxColumn% 

BoxLine$ = CHR$(186) 

PRINT BoxLine$ ; 

LOCATE BoxRow% + 1, BoxColumn% + 1 + BoxLength% 
BoxLine$ = CHR$(186) 

PRINT BoxLine$; 

LOCATE BoxRow% + 2, BoxColumn% 

BoxLine$ = CHR$(200) 

FOR it ■ 1 TO BoxLength% 

BoxLine$ = BoxLine$ + CHR$(205) 

NEXT i% 

BoxLine$ = BoxLine$ + CHR$(188) 

PRINT BoxLine$ ; 

COLOR 11, 0 
RETURN 

WriteDataToFile : 

PRINT #1, 

PRINT #1, TestID$ 
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DATABASE.BAS 


PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
PRINT #1, 
RETURN 


DataSource$ 

TestDate$ 

BumperMa t er ia 1 $ 
BumperThickness# 

Bumper StandOf f # 
PressureWallMaterial$ 
PressureWallThickness# 
Pro jectileMaterial$ 

Pro jectileDiameter# 
ImpactAngle# 

Pro jectileVelocity# 
BumperMa j or Axi s # 
BumperMinorAxi s # 
MLIHoleDiam# 
MLIMassLoss# 

Pres sWallMaj Axis# 

P r e s s Wa 1 IMi nAx i s # 


Data I npu t P r ob 1 em s 
GOTO AddDatal 
RETURN 



INVRMETH.BAS 


INVRMETH.BAS 

Source code for the inverse R method damage prediction program 

This program predicts hypervelocity impact damage using the l/R^N-l) 
interpolation/extrapolation scheme described in the MLIBlast Manual . 
Dimensioned for 100 data points and 50 active degrees of freedom. 


DIM ActiveDOF% ( 1 TO 50) 

DIM TestID$ ( 1 TO 100) 

DIM DataSource$ ( 1 TO 100) 

DIM TestDate$ ( 1 TO 100) 

DIM MaterialDataTypes$ ( 1 TO 11) 

'Dimensioned for 10 material property attributes for each material. 

DIM BumperMaterial# ( 1 TO 100, 1 TO 10) 

'Store the average of the material property attributes for the bumper. 
DIM BumperMaterialAve# ( 1 TO 10) 

DIM BumperThickness# ( 1 TO 100) 

DIM BumperStandOf f#( 1 TO 100) 

DIM PressureWallMaterial# ( 1 TO 100, 1 TO 10) 

DIM PressureWallMaterialAve# ( 1 TO 10) 

DIM PressureWallThickness# ( 1 TO 100) 

DIM Pro jectileMaterial# ( 1 TO 100, 1 TO 10) 

DIM Pro jectileMaterialAve# ( 1 TO 10) 

DIM Pro jectileDiameter# ( 1 TO 100) 

DIM ImpactAngle# ( 1 TO 100) 

DIM Pro jectileVelocity# ( 1 TO 100) 

DIM BumperMa jorAxis# ( 1 TO 100) 

DIM BumperMinorAxis# ( 1 TO 100) 

DIM MLIHoleDiam# ( 1 TO 100) 

DIM MLIMassLoss# ( 1 TO 100) 

DIM PressWallMa jAxis# ( 1 TO 100) 

DIM PressWallMinAxis#( 1 TO 100) 

DIM Material# (1 TO 10) 

DIM PredictBumpMat# ( 1 TO 10) 

DIM PredictPressWallMat# ( 1 TO 10) 

DIM PredictPro jMat#( 1 TO 10) 

DIM Xcalc#( 1 TO 10) 

DIM Fcalc#( 1 TO 10) 

DIM Weight# (1 TO 10) 

DIM a#( 1 TO 5, 1 TO 5) 

DIM B# ( 1 TO 5) 

DIM c# ( 1 TO 5) 

' Set up s green. 

COLOR 9, 0 
CLS 

LOCATE 1, 1 
PRINT CHR$ ( 201 ) ; 

FOR i% ■ 1 TO 78 

IF i% <> 39 THEN PRINT CHR$(205); 

IF i% = 39 THEN PRINT CHR$(203); 

NEXT i% 

PRINT CHR$ (187) 

LOCATE 2, Is PRINT CHR$(186) 

LOCATE 2, 40: PRINT CHR$(186) 

LOCATE 2, 80: PRINT CHR$(186) 

LOCATE 3, 1 138 

PRINT CHR$ ( 200 ) ; 



INVRMETH.BAS 


FOR 1 % = 1 TO 78 

IF i% <> 39 THEN PRINT CHR$(205); 

IF i% = 39 THEN PRINT CHR$(202); 

NEXT i% 

PRINT CHR$ ( 188 ) 

ON ERROR GOTO TestDataFileError 'Trap input test data file name errors. 

LOCATE 2, 3 
row% = 2 
col% = 24 
COLOR 11, 0 

PRINT "Test Data File Name? " ; 

COLOR 12, 0 

INPUT MLITestDataFile$ 

OPEN MLITestDataFile$ FOR INPUT AS #1 

ON ERROR GOTO MaterialDataFileError 'Trap input material data file name errors. 

LOCATE 2, 42 

row% = 2 

col% = 67 

COLOR 11, 0 

PRINT "Material Data File Name? 

COLOR 12, 0 

INPUT "", MaterialDataFile$ 

OPEN MaterialDataFile$ FOR INPUT AS #2 
ON ERROR GOTO 0 

' Count how many data values are given for each material type. 

' The user can include 10 properties on each material type. 

NumMatDat% = 0 'Stores the number of data values for each material type. 
NumMatDatl : 

INPUT #2, MaterialDataTypes$ (NumMatDat% + 1) 

IF MaterialDataTypes$ (NumMatDat% + 1) = GOTO NumMatDat2: 

NumMatDat% = NumMatDat% + 1 
GOTO NumMatDatl 
NumMatDat2 : 

CLOSE #2 

NumData% = 0 'Number of data records in the database file. 

' Intialize all the average values variables to zero. 

BumpThkAve# = 0 
BumpStandOf f Ave# » 0 
PressWallThkAve# = 0 
ProjDiaAve# = 0 
ImpAngAve# = 0 
ProjVelAve# = 0 
BumpMaj Axis Ave# = 0 
BumpMinAxisAve# = 0 
MLIHoleDiamAve# = 0 
MLIMassLossAve# = 0 
PressWallMa jAxisAve# = 0 
PressWallMinAxisAve# = 0 
'Modify top line of boarder. 

COLOR 9, 0 
LOCATE 3, 1 
PRINT CHR$( 204 ) 

LOCATE 3, 80 
PRINT CHR$ ( 185 ) 

VIEW PRINT 4 TO 12 

' Set row% to -1 to trap material name not found type of error in test data fil 
row% = -1 

DO WHILE NOT EOF(l) 

NumData% = NumData% +1 139 

INPUT #1, dummy $ 
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INPUT #1, TestID$(NumData%) 

INPUT #1, DataSource$(NumData%) 

INPUT #1, TestDate$(NumData% ) 

COLOR 9, 0 
PRINT CHR$( 186 ) ; 

COLOR 11, 0 
PRINT " No . s " ? 

PRINT NuxnData% ; 

PRINT " ID: 

PRINT TestID$ ( NumDa ta% ) ; 

PRINT " Sources " ; 

PRINT DataSource$ (NumData% ) ; 

PRINT " Date: 

PRINT TestDate$ (NumData% ) ; 

LOCATE CSRLIN, 80 
COLOR 9, 0 
PRINT CHR$ ( 186 ) 

' Input bumper material properties. 

INPUT #1, MaterialID$ 

GOSUB GetMaterialProp 'Read in material properties associated with material n 
FOR i% = 1 TO NumMatDat% 

BumperMaterial# (NumData% , i%) = Material# (i%) 

NEXT i% 

INPUT #1, BumperThickness# (NumData% ) 

BumpThkAve# - BumpThkAve# + BumperThickness# (NumData%) 

INPUT #1, BumperStandOf f # (NumData% ) 

BumpStandOf fAve# = BumpStandOf f Ave# + BumperStandOf f#(NumData% ) 

' Input pressure wall material properties. 

INPUT #1, MaterialID$ 

GOSUB GetMaterialProp 
FOR it - 1 TO NumMatDat% 

Pres sureWa llMa ter ial# (NumData% , i%) '* Material# (i%) 

NEXT i% 

INPUT #1, PressureWallThickness#( NumDa ta%) 

PressWallThkAve# = PressWallThkAve# + PressureWallThickness# (NumData% ) 

' Input projectile material properties. 

INPUT #1, Material ID$ 

GOSUB GetMaterialProp 
FOR it = 1 TO NumMatDat% 

ProjectileMaterial#(NumData%, i%) ■ Material# (i%) 

NEXT i% 

INPUT #1, Pro jectileDiameter# ( NumData% ) 

ProjDiaAve# = ProjDiaAve# + Pro jectileDiameter# (NumData% ) 

INPUT #1, ImpactAngle# (NumData%) 

ImpAngAve# » ImpAngAve# + ImpactAngle# ( NumDa ta% ) 

INPUT #1, ProjectileVelocity#(NumData%) 

ProjVelAve# * ProjVelAve# + Pro jectileVelocity#(NumData%) 

INPUT #1, BumperMa jorAxis# (NumData% ) 

BumpMa jAxisAve# = BumpMa jAxisAve# + BumperMa jorAxis# (NumData% ) 

INPUT #1, BumperMinorAxis#(NumData%) 

BumpMinAxi s Ave# = BumpMinAxisAve# + BumperMinorAxis#(NumData%) 

INPUT #1, MLIHoleDiam# ( NumData% ) 

MLIHoleDiamAve# * MLIHoleDiamAve# + MLIHoleDiam#(NumData% ) 

INPUT #1, MLIMassLoss# ( NumData% ) 

MLIMassLossAve# = MLIMassLossAve# + MLIMassLoss# (NumData%) 

INPUT #1, PressWallMa jAxis#(NumData% ) 

PressWallMa jAxisAve# = PressWallMa jAxisAve# + PressWallMa jAxis# (NumData% ) 
INPUT #1, Pres sWa 1 lMinAxi s # ( NumDa ta % ) 

PressWallMinAxisAve# = PressWallMinAxisAve# + PressWallMinAxis#( NumDa ta%) 
INPUT #1, dummy $ 


140 



invrmeth.bas 


LOOP 

VIEW PRINT 

'build box for averages 
LOCATE 12, 2 
FOR i% = 1 TO 78 
PRINT CHR$ ( 205 ) ; 

NEXT i% 

LOCATE 12, Is PRINT CHR$(204) 

LOCATE 12, 40: PRINT CHR$(203) 

LOCATE 12, 80: PRINT CHR$(185) 

FOR i% = 13 TO 23 

LOCATE i%, 1: PRINT CHR$(186) 

LOCATE i%, 40: PRINT CHR$(186) 

LOCATE i%, 80: PRINT CHR$(186) 

NEXT i% 

LOCATE 24, I: PRINT CHR$(200); 

FOR i% = 1 TO 78 

IF i% <> 39 THEN PRINT CHR$(205); 

IF i% = 39 THEN PRINT CHR$(202); 

NEXT i% 

PRINT CHR$ ( 188 ) 

' Calculate and print out database average values . 
BumpThkAve# = BumpThkAve# / NumData% 

LOCATE 13, 3 
COLOR 11, 0 

PRINT "Ave . Bumper Thk ( in ) : " ; 

COLOR 15, 0 

PRINT USING "##.### BumpThkAve# ; 

BumpS tandOff Ave# = BumpS tandOff Ave# / NumData% 
LOCATE 13, 42 
COLOR 11, 0 

PRINT "Ave. Bump. Stand Off (in):"; 

COLOR 15, 0 

PRINT USING "##.### BumpStandOf f Ave# ; 

PressWallThkAve# = PressWallThkAve# / NumData% 
LOCATE 15, 3 
COLOR 11, 0 

PRINT "Ave. Pres Wall Thk (in):"; 

COLOR 15, 0 

PRINT USING "## . ###'"''*'' ” ; PressWallThkAve#; 
ProjDiaAve# = ProjDiaAve# / NumData% 

LOCATE 15, 42 
COLOR 11, 0 

PRINT "Ave. Proj. Dia. (in):"; 

COLOR 15, 0 

PRINT USING "##.### "; ProjDiaAve#; 

ImpAngAve# = ImpAngAve# / NumData% 

LOCATE 17, 3 
COLOR 11, 0 

PRINT "Ave. Impact Angle (deg):"; 

COLOR 15, 0 

PRINT USING "##.### "; ImpAngAve#; 

ProjVelAve# = ProjVelAve# / NumData% 

LOCATE 17, 42 
COLOR 11, 0 

PRINT "Ave. Proj. Vel . (km/sec):"; 

COLOR 15, 0 

PRINT USING "##.### / ' / ' AA " ; ProjVelAve#; 

BumpMaj Axis Ave# = BumpMaj Axis Ave# / NumData% 

LOCATE 19, 3 
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COLOR 11, 0 

PRINT "Ave. Maj. Bumper Hole (in)s"; 

COLOR 15, 0 

PRINT USING "##.### AAAA "; BumpMa jAxisAve# ; 

BumpMinAxisAve# = BumpMinAxisAve# / NumData% 

LOCATE 19, 42 
COLOR 11, 0 

PRINT "Ave. Min. Bumper Hole (in):"; 

COLOR 15, 0 

PRINT USING "##.### A BumpMinAxisAve#; 

MLIHoleDiamAve# = MLIHoleDiamAve# / NumData% 

LOCATE 21, 3 
COLOR 11, 0 

PRINT "Ave. MLI Hole Diam. (in):"; 

COLOR 15, 0 

PRINT USING "##.### MLIHoleDiamAve#; 

MLIMassLossAve# = MLIMassLossAve# / NumData% 

LOCATE 21, 42 
COLOR 11, 0 

PRINT "Ave. MLI Mass Loss (grams):"; 

COLOR 15, 0 

PRINT USING "##.### A "; MLIMassLossAve#; 

Pres sWallMa jAxisAve# = PressWallMa jAxisAve# / NumData% 

LOCATE 23, 3 
COLOR 11, 0 

PRINT "Ave. Maj. P.Wall Hole (in):"; 

COLOR 15, 0 

PRINT USING "##.###'' AAA " ; PressWallMa jAxisAve#; 

PressWallMinAxisAve# = PressWallMinAxisAve# / NumData% 

LOCATE 23, 42 
COLOR 11, 0 

PRINT "Ave. Min. P.Wall Hole (in):"; 

COLOR 15, 0 

PRINT USING "##.### AAAA " ; PressWallMinAxisAve#; 

LOCATE 25, 28 
COLOR 12, 0 

PRINT "press any key to continue"; 

DO 

LOOP WHILE INKEY$ = 'Press any key to continue 

f 

• Normalize data with respect to calculated averages and. 

' determine which DOF are active. ActiveDOF%(i%) = 0 if i-th 
' DOF is not active and = 1 if is active (ie - changes in the database). 


FOR i% = 1 TO 50 

ActiveDOF% ( i% ) = 0 
NEXT i% 

FOR it = 1 TO NumData% 

BumperThickness# ( i% ) = BumperThickness#(i%) / BumpThkAve# 

IF BumperThickness# (i%) < .999 OR BumperThickness# (i%) > 1.001 THEN ActiveDOF 
BumperStandOf f # ( i% ) = BumperStandOf f #(i% ) / BumpStandOf fAve# 

IF BumperStandOf f# ( i% ) < .999 OR BumperStandOf f#(i%) > 1.001 THEN ActiveDOF%( 
PressureWallThickness# ( i% ) = PressureWallThickness# (i% ) / PressWallThkAve# 

IF PressureWallThickness# (i%) < .999 OR PressureWallThickness# ( i% ) > 1.001 TH 
Pro jectileDiameter# ( i% ) = Pro jectileDiameter#(i%) / ProjDiaAve# 

IF Pro jectileDiameter# ( i% ) < .999 OR ProjectileDiameter#(i%) > 1.001 THEN Act 
IF ImpAngAve# <> 0 THEN ImpactAngle# ( i% ) = ImpactAngle# (i% ) / ImpAngAve# 

IF ImpactAngle# (i%) < .999 OR ImpactAngle# (i%) > 1.001 THEN ActiveDOF% ( 5 ) - 1 
ProjectileVelocity#(i%) = Pro jectileVelocity#(i%) / ProjVelAve# 

IF Pro jectileVelocity# ( i% ) < .999 OR ProjectileVelocity#(i%) > 1.001 THEN Act 
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NEXT i% 

' Normalize the material data 
FOR i% = 1 TO NumMatDat% 

FOR j% - 1 TO NumData% 

BumperMaterialAve# ( i% ) = BumperMaterialAve# (i%) + BumperMaterial# ( j% , i%) 
PressureWallMaterialAve#(i%) = PressureWallMaterialAve#(i%) + PressureWall 
Pro jectileMaterialAve# ( i% ) = ProjectileMaterialAve#(i%) + Pro jectileMateri 
NEXT j% 

NEXT i% 

FOR i% = 1 TO NumMatDat% 

BumperMaterialAve# (i%) = BumperMaterialAve# (i%) / NumData% 
PressureWallMaterialAve# ( i% ) = PressureWallMaterialAve# ( i% ) / NumData% 

Pro jectileMaterialAve# ( i% ) = Pro jectileMaterialAve# (i%) / NumData% 


NEXT i% 

' Display material property 
CLS 

LOCATE 2, 1 
COLOR 11, 0 

information. 



PRINT " 

MATERIAL 

BUMPER 

PRESSURE 

PROJECT 

PRINT " 

PROPERTY 

MATERIAL 

WALL 

MATERI 

PRINT " 

FOR i% = 1 TO 
LOCATE 2 * 

NumMatDat% 
i% + 4, 3 

AVE. 

MATERIAL AVE. 

AVE 


COLOR 11, 0 

PRINT MaterialDataTypes$ ( i% ) 

LOCATE 2 * i% + 4, 33 
COLOR 15, 0 

PRINT USING "##.###. BumperMaterialAve# (i%) 

LOCATE 2 * i% + 4, 50 

PRINT USING "##.### PressureWallMaterialAve# ( i% ) 

LOCATE 2 * i% + 4, 67 

PRINT USING "##.### Pro jectileMaterialAve# ( i% ) ; 

NEXT i% 

BottomRow% = CSRLIN 
' Draw Box Around The Data 
FOR i% = 2 TO 79 
COLOR 9, 0 

LOCATE 1, i% : PRINT CHR$(205) 

LOCATE 5, i% : PRINT CHR$(205) 

LOCATE BottomRow% + 1, i%: PRINT CHR$(205); 

NEXT i% 

FOR i% = 2 TO BottomRow% 

LOCATE i%, 1: PRINT CHR$(186) 

LOCATE i% , 29: PRINT CHR$(186) 

LOCATE i% , 46: PRINT CHR$(179) 

LOCATE i% , 63: PRINT CHR$(179) 

LOCATE i% , 80: PRINT CHR$(186) 

NEXT i% 

LOCATE 1, 1: PRINT CHR$(201) 

LOCATE 1, 29: PRINT CHR$(203) 

LOCATE 1, 46: PRINT CHR$(209) 

LOCATE 1, 63: PRINT CHR$(209) 

LOCATE 1, 80: PRINT CHR$(187) 

LOCATE 5, 1: PRINT CHR$(204) 

LOCATE 5, 29: PRINT CHR$(206) 

LOCATE 5, 46: PRINT CHR$(216) 

LOCATE 5, 63: PRINT CHR$(216) 

LOCATE 5, 80: PRINT CHR$(185) 

LOCATE BottomRow% +1, 1: PRINT CHR$(200); 

LOCATE BottomRow% + 1, 29: PRINT CHR$(202); 
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LOCATE BottomRow% + 1, 46: PRINT CHR$(207); 

LOCATE BottomRow% + 1, 63: PRINT CHR$(207); 

LOCATE BottomRow% + 1, 80: PRINT CHR$(188); 

LOCATE 25, 28 
COLOR 12, 0 

PRINT "press any key to continue"; 

DO 

LOOP WHILE INKEY$ = "" 'Press any key to continue 

' Check which material properties are active in the database. 

FOR i% = 1 TO NumMatDat% 

FOR 14 = 1 TO NumData% 

BumperMaterial# ( j % , i%) = BumperMaterial#* j%, i%) / BumperMatenalAve#(i%) 
IF BumperMaterial# ( j % , i%) < .999 OR BumperMaterial# ( j% , i%) > 1.001 THEN 
PressureWallMaterial#* j%, i%) = PressureWallMaterial#* j%, i%) / PressureWa 
IF PressureWallMaterial# (j%, i%) < .999 OR PressureWallMaterial# ( j%, i%) > 
Pro jectileMaterial# ( j% , i%) = Pro jectileMaterial#( j%, i%) / Pro jectileMate 
IF Pro jectileMaterial# ( j%, i%) < .999 OR Pro jectileMaterial# ( j%, i%) > 1.0 
NEXT j% 

NEXT i% 

NumActiveDOF% =0 _ . 

FOR i% = 1 TO 6 + 3 * NumMatDat% 'Count how many degrees of freedom are active. 

NumActiveDOF% = NumActiveDOF% + ActiveDOF%(i%) 

NEXT i% 

' Set the default values. These can be changed. 

BumperMaterial ID$ “ "6061-T6" ... ...... .. 

PredictBumpThick# = . 04 / BumpThkAve# 'Default values should be divided by the 
PredictBumpStandOf f # = 4 / BumpS tandOffAve# 

PressWallMaterialID$ = "2219-T87" 

PredictPressWallThick# = .125 / PressWallThkAve# 

ProjectileMaterialID$ = "1100" 

PredictPro j Di a # = .313 / ProjDiaAve# 

IF ImpAngAve# <> 0 THEN PredictlmpAngle# = 45 / ImpAngAve# 

PredictPro jVel# = 5.3# / ProjVelAve# 

t 

' Prompt the user for the prediction required. 


PredictValue : 

SCREEN 0 
COLOR 15, 1 
CLS 

LOCATE 1, 3 

PRINT "Please Enter The "; 

COLOR 12, 1 
PRINT "Number 
COLOR 15, 1 

PRINT "Associated With the Desired Action:" 
LOCATE 5, 3 
COLOR 12, 1 
PRINT " 1 . " ; 

COLOR 15, 1 

PRINT "Predict Bumper Hole Major Axis" 
LOCATE 7, 3 
COLOR 12, 1 
PRINT " 2 . " ; 

COLOR 15, 1 

PRINT "Predict Bumper Hole Minor Axis" 
LOCATE 9, 3 
COLOR 12, 1 
PRINT "3. "; 


COLOR 15, 1 
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PRINT "Predict MLI Hole Diameter" 

LOCATE 11, 3 
COLOR 12, 1 
PRINT "4. "? 

COLOR 15, 1 

PRINT "Predict MLI Mass Loss" 

LOCATE 13, 3 
COLOR 12, 1 
PRINT " 5 . " ? 

COLOR 15, 1 

PRINT "Predict Pressure Wall Hole Major Axis" 

LOCATE 15, 3 
COLOR 12, 1 
PRINT " 6 . " ; 

COLOR 15, 1 

PRINT "Predict Pressure Wall Hole Minor Axis" 

SelectPredictionType : 

LOCATE 19, 1 

COLOR 12, 1 

INPUT PredictionType% 

IF PredictionType% < 1 OR PredictionType% > 6 THEN 
LOCATE 21, 1 

PRINT "Please re-enter choice I" 

LOCATE 19, 1 
PRINT " 

GOTO SelectPredictionType 
END IF 

f 

' Prompt user for data associated with prediction. 

l 

COLOR 11, 0 
CLS 

LOCATE 1, 1 

PRINT "ENTER DATA FOR DESIRED PREDICTION! " 

COLOR 15, 0 
LOCATE 3, 1 

PRINT "[default values shown in square brackets]" 

COLOR 10, 0 
LOCATE 5, 1 

PRINT "(magnitude relative to database average shown in round brackets)" 

LOCATE 8, 1 
COLOR 11, 0 

PRINT "Bumper Material!"; 

COLOR 15, 0 

PRINT " [ " ; BumperMaterialID$ ; " ] " 

LOCATE 8, 40 
row% ■ 8: col% = 40 
COLOR 12, 0 
INPUT " dummy$ 

IF dummy$ “ " " THEN 

Material ID$ = BumperMaterialID$ 

ELSE 

BumperMaterialID$ = dummy$ 

Material ID$ = dummy $ 

END IF 

Material ID$ = UCASE$ (Material ID$ ) 

GOSUB GetMaterialProp 
BumperMaterialID$ = Material ID$ 

MatAve# = 0 'Used to keep track of average material properties. 

DOFWaming% = 0 'A parameter used to warn user that a prediction is requested 
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'for a DOF that does not vary in the database. 
FOR it a 1 TO NumMatDat% 

Predi C tBui n pMa t# ( i %) = Material# ( i% ) / BumperMaterialAve# ( i% ) 
IF (Predi c tBumpM a t#(i%) < .999 OR PredictBumpMat#(i%) > 1.00 
Mai- lira* = MatAve# + PredictBumpMat# ( i% ) 


MatAve# 
NEXT i% 


001) AND ActiveDOF% 


MatAve# = MatAve# / NumMatDat% 

LOCATE 8, 60 
COLOR 10, 0 
PRINT 

PRINT USING ”##.###"; MatAve#; 

PRINT " ) " 

IF DOFWarning% = 1 THEN 
LOCATE 9, 1 
COLOR 12, 0 

Warning-Your test data does not support the above input data!” 

hNU IF 


DOFWarning% = 0 
LOCATE 10, 1 

PredictBuinpThick# = PredictBumpThick# * BumpThkAve# 

COLOR 11, 0 

PRINT "Bumper Thickness (in):"; 

COLOR 15, 0 
PRINT " [ " ; 

PRINT USING "##.###'■''*''"; PredictBumpThick#; 

PRINT "]" 

LOCATE 10, 40 
COLOR 12, 0 
INPUT dummy# 

IF dummy# <> 0 THEN PredictBumpThick# = dummy# 

PredictBumpThick# = PredictBumpThick# / BumpThkAve# 

LOCATE e i0 Ct 60 mPThiCk# < ’" 9 ° R PredictBum P Thick * > 1-001) AND ActiveDOF% ( 1 ) = 0 

COLOR lo/o 
PRINT "(" ; 

PRINT USING ”##.###"; PredictBumpThick#; 

PRINT ") M 

IF DOFWarning% = 1 THEN 
LOCATE 11, 1 
COLOR 12, 0 

Warnin 9“ Your test data does not support the above input data I" 


IF Sralngl^! >- 5 THEN ' Include tha ®«ects o£ the pressure wall here. 
LOCATE 12, 1 

PredictBumpS tandOf f # = PredictBumpStandOf f # * BumpStandOf fAve# 

COLOR 11,0 

PRINT "Bumper Stand-off (in):"; 

COLOR 15, 0 
PRINT " [ " ; 

PRINT USING "##.###"■""''" ; PredictBumpStandOf f# ; 

PRINT "]" 

COLOR 12, 0 
LOCATE 12, 40 
INPUT dummy# 

IF dummy# <> 0 THEN PredictBumpStandOf f# = dummy# 

PredictBumpStandOf f# = PredictBumpStandOf f# / BumpStandOf fAve# 

ro^lmo e fo Ct f!f npStanCi0ff * < '" 9 0R Predic tBumpStandOff# > 1.001) AND ActiveDO 
LOCATE 12 f 60 

COLOR 10, 0 146 



PRINT " ( " ; 

PRINT USING "##.###" ; PredictBumpStandOf f #? 

PRINT " ) " 

IF DOFWaming% = 1 THEN 
LOCATE 13, 1 
COLOR 12, 0 

PRINT "Warning-Your test data does not support the above input data l " 

END IF 

DOFWarning% = 0 
LOCATE 14, 1 
COLOR 11, 0 

PRINT "Pressure Wall Materials"; 

COLOR 15, 0 

PRINT " [ " ; PressWallMaterialID$ ; " ] " 

LOCATE 14, 40 
row% = 14: col% = 40 
COLOR 12, 0 
INPUT " " , dummy? 

IF dummy? = "" THEN 

MateriallD? = PressWallMateriallD? 

ELSE 

PressWallMateriallD? = dummy? 

MateriallD? = dummy? 

END IF 

MateriallD? = UCASE? (MateriallD? ) 

GOSUB GetMaterialProp 
PressWallMateriallD? = MateriallD? 

MatAve# = 0 

FOR i% = 1 TO NumMatDat% 

PredictPressWallMat# ( i% ) = Material# ( i% ) / PressureWallMaterialAve# ( i% ) 

IF (PredictPressWallMat# (i%) < .999 OR PredictPressWallMat# ( i% ) > 1.001) A 
MatAve# = MatAve# + PredictPressWallMat# (i%) 

NEXT i% 

MatAve# = MatAve# / NumMatDat% 

LOCATE 14, 60 
COLOR 10, 0 
PRINT " ( " ; 

PRINT USING "##.###"; MatAve#; 

PRINT ")" 

IF DOFWaming% = 1 THEN 
LOCATE 15, 1 
COLOR 12, 0 

PRINT "Warning-Your test data does not support the above input data i " 

END IF 

DOFWarning% = 0 
LOCATE 16, 1 

PredictPressWallThick# = PredictPressWallThick# * PressWallThkAve# 

COLOR 11, 0 

PRINT "Press. Wall Thick, (in):"; 

COLOR 15, 0 
PRINT "["; 

PRINT USING "##.###''~''' s "; PredictPressWallThick#; 

PRINT "]" 

LOCATE 16, 40 
COLOR 12, 0 
INPUT "", dummy# 

IF dummy# <> 0 THEN PredictPressWallThick# = dummy# 

PredictPressWallThick# = PredictPressWallThick# / PressWallThkAve# 

IF (PredictPressWallThick# < .999 OR PredictPressWallThick# > 1.001) AND Acti 
LOCATE 16, 60 1/7 
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COLOR 10, 0 
PRINT 

PRINT USING "##.###"; PredictPressWallThick# ; 

PRINT " ) H 

IF DOFWaming% = 1 THEN 
LOCATE 17, 1 
COLOR 12, 0 

PRINT "Warning-Your test data does not support the above input data l" 

END IF 
END IF 

DOFWarning% ” 0 
LOCATE 18, 1 
COLOR 11, 0 

PRINT "Projectile Material:"; 

COLOR 15, 0 

PRINT " [ " ; Pro jectileMaterialID$ ; " ] " 

LOCATE 18, 40 
row% * 18: col% = 40 
COLOR 12, 0 
INPUT " ", dummy$ 

IF dummy$ = "" THEN 

MaterialID$ ■ Pro jectileMaterialID$ 

ELSE 

Pro jectileMaterialID$ = dummy$ 

MaterialID$ • dumray$ 

END IF 

MaterialID$ = UCASE$ (MaterialID$ ) 

GOSUB GetMaterialProp 
ProjectileMaterialID$ = MaterialID$ 

MatAve# = 0 

FOR i% =» 1 TO NumMatDat% 

PredictPro jMat#(i% ) = Material# (i%) / ProjectileMaterialAve#(i%) 

IF ( PredictPro jMat#(i%) < .999 OR PredictPro jMat#(i%) > 1.001) AND ActiveDOF 1 
MatAve# = MatAve# + PredictPro jMat# ( i% ) 

NEXT i% 

MatAve# = MatAve# / NumMatDat% 

LOCATE 18, 60 
COLOR 10, 0 
PRINT "("; 

PRINT USING "##.###"; MatAve#; 

PRINT ")" 

IF DOFWarning% = 1 THEN 
LOCATE 19, 1 
COLOR 12, 0 

PRINT "Waming-Your test data does not support the above input data ! " 

END IF 

DOFWarning% = 0 
LOCATE 20, 1 

PredictPro jDia# = PredictPro jDia# * ProjDiaAve# 

COLOR 11, 0 

PRINT "Projectile Diameter (in):"; 

COLOR 15, 0 
PRINT "["; 

PRINT USING ”##.###^ A/k "; PredictPro jDia#; 

PRINT "]" 

LOCATE 20, 40 
COLOR 12, 0 
INPUT "", dummy# 

IF dummy# <> 0 THEN PredictPro jDia# - dummy# 

PredictPro jDia# » PredictPro jDia# / ProjDiaAve# 
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IF ( PredictPro jDia# < .999 OR PredictPro jDia# > 1.001) AND ActiveDOF% ( 4 ) 3J THE 
LOCATE 20, 60 
COLOR 10, 0 
PRINT 

PRINT USING "##.###"; PredictPro jDia#; 

PRINT " ) " 

IF DOFWarning% = 1 THEN 
LOCATE 21, 1 
COLOR 12, 0 

PRINT "Warning-Your test data does not support the above input data ! " 

END IF 

DOFWarning% = 0 
LOCATE 22, 1 
COLOR 11, 0 

PRINT "Impact Angle (degrees):"; 

LOCATE 22, 40 
COLOR 12, 0 

INPUT "", Predict ImpAngle# 

IF ImpAngAve# <> 0 THEN Predict ImpAngle# = Predict ImpAngle# / ImpAngAve# 

IF ImpAngAve# = 0 AND Predict ImpAngle# > 0 THEN 
LOCATE 23, 1 

COLOR 12, 0 . 

PRINT "Warning-Your test data does not support the above input data l 

END IF 

IF ImpAngAve# <> 0 THEN 

IF (Predict ImpAngle# < .999 OR PredictlmpAngle# > 1.001) AND ActiveDOF% ( 5 ) = 
LOCATE 22, 60 
COLOR 10, 0 
PRINT " ( " ; 

PRINT USING "##.###"; PredictlmpAngle#; 

PRINT " )" 

IF DOFWaming% = 1 THEN 
LOCATE 23, 1 
COLOR 12, 0 

PRINT "Warning-Your test data does not support the above input data!" 

END IF 
END IF 

DOFWarning% = 0 
LOCATE 24 , 1 

PredictPro jVel# = PredictPro jVel# * ProjVelAve# 

COLOR 11, 0 

PRINT "Proj. Vel . (km/sec):"; 

COLOR 15, 0 
PRINT "["; 

PRINT USING ”##.###~ / ' < ' / '"; PredictPro jVel#; 

PRINT " ]"; 

LOCATE 24, 40 
COLOR 12, 0 
INPUT ; "", dummy# 

IF dummy# <> 0 THEN PredictPro jVel# = dummy# 

PredictPro j Vel# = PredictPro j Vel# / ProjVelAve# 

IF (PredictPro jVel# < .999 OR PredictPro j Vel# > 1.001) AND ActiveDOF% ( 6 ) = 0 THE 
LOCATE 24, 60 
COLOR 10, 0 
PRINT "("; 

PRINT USING "##.###"; PredictPro j Vel#; 

PRINT ")"; 

IF DOFWarning% = 1 THEN 

LOCATE 25, 1 1^9 

COLOR 12, 0 
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PRINT "Waming-Your test data does not support the above input data ! " ; 

END IF 

LOCATE 25, 28 
COLOR 12, 0 

PRINT "press any key to continue"; 

DO 

LOOP WHILE INKEY$ ■ "" 'Press any key to continue 

r 

' Find distance of user input point from origin. 

f 

PredictVectMag# = 0 
FOR i% = 1 TO NumMatDat% 

PredictVectMag# = PredictVectMag# + PredictBumpMat#(i%) ~ 2 
IF PredictionType% >= 5 THEN 

PredictVectMag# * PredictVectMag# + PredictPressWallMat#(i%) A 2 
END IF 

PredictVectMag# = PredictVectMag# + PredictProjMat#(i%) ~ 2 
NEXT i% 

PredictVectMag# = PredictVectMag# + PredictBumpThick# A 2 
IF PredictionType% >= 5 THEN 

PredictVectMag# = PredictVectMag# + PredictBumpStandOf f # A 2 
PredictVectMag# = PredictVectMag# + PredictPressWallThick# " 2 
END IF 

PredictVectMag# = PredictVectMag# + PredictPro jDia# * 2 
PredictVectMag# “ PredictVectMag# + Predict ImpAngle# A 2 
PredictVectMag# = PredictVectMag# + PredictPro jVel# ~ 2 
PredictVectMag# = PredictVectMag# A . 5 

' Determine near and far points of data points along user input prediction vect 
Far# =0 

DistanceAve# = 0 

FOR i% - 1 TO NumData% 

Distance# = 0 

FOR j% = 1 TO NumMatDat% 

Distance# * Distance# + PredictBumpMat# ( j% ) * BumperMaterial#(i%, j%) 

IF PredictionType% >= 5 THEN 

Distance# = Distance# + PredictPressWallMat#( j%) * PressureWallMaterial 
END IF 

Distance# = Distance# + PredictPro jMat#( j%) * ProjectileMaterial#(i%, j%) 
NEXT j% 

Distance# = Distance# + PredictBumpThick# * BumperThickness#(i%) 

IF PredictionType% >= 5 THEN 

Distance# = Distance# + PredictBumpStandOf f# * BumperStandOf f # ( i% ) 
Distance# = Distance# + PredictPressWallThick# * PressureWallThickness#(i% 
END IF 

Distance# = Distance# + PredictPro jDia# * Pro jectileDiameter# ( i% ) 

Distance# = Distance# + Predict ImpAngle# * ImpactAngle#(i%) 

Distance# « Distance# + PredictPro jVel# * ProjectileVelocity# ( i% ) 

Distance# = Distance# / PredictVectMag# 

IF Far# = 0 THEN 

Near# = Distance# 

Far# * Distance# 

ELSE 

IF Distance# < Near# THEN Near# = Distance# 

IF Distance# > Far# THEN Far# = Distance# 

END IF 

DistanceAve# = DistanceAve# + Distance# 

NEXT i% 

DistanceAve# - DistanceAve# / NumData% 

SCREEN 0 

COLOR 9, 0 liU 
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CLS 

LOCATE 1, 1: PRINT CHR$(201) 

FOR i% - 2 TO 79 

LOCATE 1, i% : PRINT CHR$(205) 

NEXT i% 

LOCATE 1, 80s PRINT CHR$(187) 

Relatives = 100 * PredictVectMag# / DistanceAve# 

PRINT CHR$ (186) ; 

ZOLOR 11, 0 

PRINT " Relative distance of desired prediction from origin: 

COLOR 15, 0 
PRINT Relatives 
COLOR 9, 0 

LOCATE 2, 80: PRINT CHR$(186) 

LOCATE 3, 1: PRINT CHR$(186): LOCATE 3, 80: PRINT CHR$(186) 

Relatives = 100 
PRINT CHR$ ( 186 ) ; 

COLOR 11, 0 

PRINT H Relative distance of mean data point from the origin: Relatives 

COLOR 9, 0 

LOCATE 4, 80: PRINT CHR$(186) 

LOCATE 5, 1: PRINT CHR$(1 06 ) : LOCATE 5, 80: PRINT CHR$(186) 

Relatives = 100 * Far# / DistanceAve# 

PRINT CHR$ ( 186 ) ; 

COLOR 11, 0 

PRINT " Relative dist of farthest data point from the origin: Relatives 

COLOR 9, 0 

LOCATE 6, 80: PRINT CHR$(186) 

LOCATE 7, 1: PRINT CHR$(186): LOCATE 7, 80: PRINT CHR$(186) 

Relatives = 100 * Near# / DistanceAve# 

PRINT CHR$ (186 ) ; 

COLOR 11, 0 

PRINT " Relative dist of nearest data point from the origin: Relatives 

COLOR 9, 0 

LOCATE 8, 80: PRINT CHR$(186) 

LOCATE 9, 1: PRINT CHR$(204) 

FOR i% = 2 TO 79 

LOCATE 9, i% : PRINT CHR$(205) 

NEXT i% 

LOCATE 9, 80: PRINT CHR$(185) 

' Calculate ten data points along user input vector 
FOR i% = 1 TO 10 
Normalize# = 0 

Delta# = (Far# - Near#) / 20# 

Factor# = ( ( ( i% - 1) * 2# + 1#) * Delta# + Near#) / PredictVectMag# 

Xcalc#(i% ) = Factor# * PredictVectMag# 

Fcalc#(i%) = 0 
FOR )% 3 1 TO NumData% 
r# = 0 

FOR k% * 1 TO NumMatDat% 

r# = r# + (BumperMaterial#( j%, k%) - PredictBumpMat#(k%) * Factor#) A 2 
IF PredictionType% >= 5 THEN 

r# = r# + (PressureWallMaterial#( j%, k%) - PredictPressWallMat#(k%) 
END IF 

r# = r# + (Pro jectileMaterial#( j%, k%) - PredictProjMat#(k%) * Factor#) 
NEXT k% 

r# = r# + (BumperThickness#( j%) - PredictBumpThick# * Factor#) A 2 
IF PredictionType% >* 5 THEN 

r# = r# + ( Bumpers tandOf f # ( j% ) - PredictBumpStandOf f # * Factor#) A 2 
r# = r# + (PressureWallThickness#( j%) - PredictPressWallThick# * Factor 
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END IF 

r# = r # + (ProjectileDiameter#( j%) - PredictProjDia# * Factor#) " 2 
r# * r# + ( ImpactAngle# ( j % ) - PredictlmpAngle# * Factor#) A 2 
r # = r # + (Pro jectileVelocity# ( j% ) - PredictProjVel# * Factor#) * 2 
r# - r# ^ .5 

Normalize# ■ Normalize# + 1# / r# A NumActiveDOF% 

IF PredictionType% = 1 THEN Fcalc#(i%) = Fcalc#(i%) + BumperMa jorAxis# ( j% 

IF PredictionType% = 2 THEN Fcalc#(i%) = Fcalc#(i%) + BumperMinorAxis# ( j% 

IF PredictionType% - 3 THEN Fcalc#(i%) = Fcalc#(i%) + MLIHoleDiam# ( j% ) / 

IF PredictionType% = 4 THEN Fcalc#(i%) = Fcalc#(i%) + MLIMassLoss# ( j% ) / 

IF Predict ionType% - 5 THEN Fcalc#(i%) = Fcalc#(i%) + PressWallMa jAxis#( j 

IF PredictionType% - 6 THEN Fcalc#(i%) = Fcalc#(i%) + PressWallMinAxis# ( j 

NEXT j% 

IF Normalize# <> 0 THEN Fcalc#(i%) = Fcalc#(i%) / Normalize# 

NEXT i% 

f 

' Here we fit a polynomial of order Order% through interpolated data points . 

Order% ■ 5 

NumGoodData% = 0 

sumx# = 0 

sumx2# = 0 

sumx 3# =* 0 

sumx 4# = 0 

sumx 5# 3 0 

sumx 6# = 0 

sumx 7# = 0 

sumx8# = 0 

ybar# = 0 

sumy# = 0 

sumyx# = 0 

sumyx2# =* 0 

sumyx 3# = 0 

sumyx4# = 0 

SumWeight# = 0 

FOR i% - 1 TO 10 

Weight# (i%) = 1# 'Here a provision is made for weighting the data. 

SumWeight# ■ SumWeight# + Weight# (i%) 

sumx# = sumx# + Xcalc#(i%) * Weight# (i%) 

sumx2# = sumx 2# + Xcalc#(i%) A 2 * Weight# (i%) 

sumx3# * sumx3# + Xcalc#(i%) * 3 * Weight# (i%) 

sumx4# - sumx4# + Xcalc#(i%) A 4 * Weight# (i%) 

sumx5# ■ sumx5# + Xcalc#(i%) A 5 * Weight# (i%) 

sumx6# = sumx6# + Xcalc#(i%) ^ 6 * Weight# (i%) 

sumx7# = sumx7# + Xcalc#(i%) A 7 * Weight# (i%) 

sumx8# ■ sumx8# + Xcalc#(i%) * 8 * Weight# (i%) 

ybar# = ybar# + Fcalc#(i%) 

sumy# = sumy# + Fcalc#(i%) * Weight# (i%) 

sumyx# * sumyx# + Fcalc#(i%) * Xcalc#(i%) * Weight# (i%) 

sumyx2# * sumyx2# + Fcalc#(i%) * Xcalc#(i%) * 2 * Weight# (i%) 

sumyx 3# = sumyx 3# + Fcalc#(i%) * Xcalc#(i%) ^ 3 * Weight# (i%) 

sumyx4# * sumyx4# + Fcalc#(i%) * Xcalc#(i%) A 4 * Weight# (i%) 

NEXT i% 

' Here the least square equation [a]{b> = {c> is set up and solved. 

' {b> are the polynomial coefficients, 

ybar# ■ ybar# / 10 
a#(l / 1) = SumWeight# 
a#(l, 2) = sumx# 
a#^!, 3) = sumx2# 
a#(l, 4) = sumx3# 
a#(l, 5) = sumx4# 
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c#(l) * sumy# 
a#(2, 1) = sumx# 
a# (2, 2) = sumx2# 
a#(2, 3) = sumx3# 
a#(2, 4) = suxnx4# 
a#(2, 5) = sumx5# 
c# ( 2 ) = sumyx# 
a#(3, 1) = sumx2# 
a#(3, 2) = sumx3# 
a# (3, 3) = sumx4# 
a#(3, 4) = sumx5# 
a#(3, 5) = sumx6# 
c#(3) = sumyx 2# 
a#(4, 1) = sumx3# 
a#(4, 2) - sumx4# 
a#(4, 3) = sumx5# 
a#(4, 4) = sumx6# 
a#(4, 5) = sumx7# 
c#(4) = sumyx 3# 
a#(5, 1) = sumx4# 
a#(5, 2) = sumx5# 
a#(5, 3) = sumx6# 
a#(5, 4) = sumx7# 
a#(5, 5) = sumx8# 
c#(5) ■ sumyx 4# 

FOR i% ■ 1 TO Order% - 1 

FOR j% = i% + 1 TO Order% 

IF a# ( j% , i% ) <= a# ( i% , i%) THEN GOTO LI 
FOR k% = i% TO Order % 
dumb# = a#(i%, k%) 
a#(i%, k%) = a#(j%, k%) 
a#( j%, k% ) = dumb# 

NEXT k% 

dumb# = c#(i%) 
c#(i%) = c#( j%) 
c#(j%) = dumb# 

Lis 

NEXT j% 

FOR j% = i% + 1 TO Order% 

Factor# = a#( j%, i%) / a#(i%, i%) 

FOR k% = i% TO Order% 

a#( j%, k% ) = a#(j%, k% ) - a#(i%, k%) * Factor# 
NEXT k% 

c#(j%) = c#(j%) - c#(i%) * Factor# 

NEXT j% 

NEXT i% 

FOR i% = Order % TO 1 STEP -1 
sum# “ 0 

IF i% = Order% THEN 

B#(i%) = (c#(i%) - sum#) / a#(i%, i%) 

GOTO L2 
END IF 

FOR j% = i% + 1 TO Order% 

sum# = sum# + a#(i%, j%) * B#(j%) 

NEXT j% 

B# ( i % ) = (c#(i%) - sum#) / a#(i%, i%) 

L2s 
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FOR i% = 1 TO 10 

fest# =* B# ( 1 ) + B#( 2 ) * Xcalc# ( i% ) + B#(3) * Xcalc#(i%) * 2 + B#(4) * Xcalc#( 
sumerr# =* sumerr# + (Fcalc#(i%) - fest#) A 2 * Weight# (i%) 
sumdif# ■ sumdif# + (Fcalc#(i%) - ybar#) A 2 * Weight#(i%) 

NEXT i% 

IF sumdif# <> 0 THEN CoefDet# = 1 - (sumerr# / sumdif#) 

IF sumdif# = 0 THEN CoefDet# = 1 
LOCATE 10, 1 
COLOR 9, 0 
PRINT CHR$ ( 186 ) ; 

COLOR 11, 0 

PRINT " LEAST SQUARES POLYNOMIAL FIT THROUGH TEN INTERPOLATED DATA POINTS:" 

COLOR 9, 0 

LOCATE 10, 80: PRINT CHR$(186) 

COLOR 11, 0 

LOCATE 11, 3: PRINT " ( " ; NumData%; " Data Records From 
COLOR 15, 0 

PRINT MLITestDataFile$; 

COLOR 11, 0 

PRINT " Database Were Used For Interpolation)" 

COLOR 9, 0 

LOCATE 11, 1: PRINT CHR$(186): LOCATE 11, 80: PRINT CHR$(186) 

PRINT CHR$ ( 186 ) ; 

COLOR 11, 0 

PRINT " constant term " ; 

PRINT USING "##.### "; B#(l) 

COLOR 9, 0 

LOCATE 12, 80: PRINT CHR$(186) 

PRINT CHR$( 186 ) ; 

COLOR 11, 0 

PRINT " x coefficient 

PRINT USING "##.### AAAA "? B# ( 2 ) 

COLOR 9, 0 

LOCATE 13, 80: PRINT CHR$(186) 

PRINT CHR$ ( 186 ) ; 

COLOR 11, 0 

PRINT " x A 2 coefficient "; 

PRINT USING "##.### B#(3) 

COLOR 9, 0 

LOCATE 14, 80: PRINT CHR$(186) 

PRINT CHR$( 186 ) ; 

COLOR 11, 0 

PRINT " x A 3 coefficient "; 

PRINT USING "##.### "; B#(4) 

COLOR 9, 0 

LOCATE 15, 80: PRINT CHR$(186) 

PRINT CHR$ (186); 

COLOR 11, 0 

PRINT " x A 4 coefficient 

PRINT USING "##.### B# ( 5 ) 

COLOR 9, 0 

LOCATE 16, 80: PRINT CHR$(186) 

PRINT CHR$ (186 ) ; 

COLOR 11, 0 

PRINT ■ POLYNOMIAL coefficient of determination (R A 2) *> 

COLOR 15, 0 

PRINT USING "## •### / ' AAA " ; CoefDet# 

COLOR 9, 0 

LOCATE 17, 80: PRINT CHR$(186) 

LOCATE 18, 1: PRINT CHR$(204) 


154 



INVRMETH.BAS 

FOR i% = 2 TO 79 

LOCATE 18, i% : PRINT CHR$(205) 

NEXT i% 

LOCATE 18, 80: PRINT CHR$(185) 

LOCATE 19, 1: PRINT CHR$(186): LOCATE 19, 80: PRINT CHR$(186) 

PREDICTION# = B# ( 1 ) + B# ( 2 ) * PredictVectMag# + B#(3) * PredictVectMag# A 2 + B# 
IF PREDICTION# < 0 THEN PREDICTION# = 0 
PRINT CHR$ ( 186 ) ; 

COLOR 11, 0 

IF Predict ionType% = 1 THEN PRINT " Predicted Bumper Hole Major Axis (in) ="; 

IF PredictionType% = 2 THEN PRINT " Predicted Bumper Hole Minor Axis (in) ="; 

IF PredictionType% = 3 THEN PRINT " Predicted MLI Hole Diameter (in) ="; 

IF PredictionType% = 4 THEN PRINT " Predicted MLI Mass Loss (grams) = "; 

IF PredictionType% = 5 THEN PRINT " Predicted Pressure Wall Hole Major Axis (in) 

IF PredictionType% = 6 THEN PRINT " Predicted Pressure Wall Hole Minor Axis (in) 

COLOR 15, 0 

PRINT USING "##.### PREDICTION# 

COLOR 9, 0 

LOCATE 20, 80: PRINT CHR$(186) 

LOCATE 21, 1: PRINT CHR$(186): LOCATE 21, 80: PRINT CHR$(186) 

LOCATE 22, 1: PRINT CHR$(200) 

FOR i% = 2 TO 79 

LOCATE 22, i%: PRINT CHR$(205) 

NEXT i% 

LOCATE 22, 80: PRINT CHR$(188) 

LOCATE 24, 15 
COLOR 9, 0 

PRINT "Do you wish to see a plot of the results (y/n)? 

COLOR 12, 0 
INPUT " ", Answer? 

COLOR 15, 0 

Answer? = LCASE? (Answer? ) 

IF Answer? = "y" THEN GOSUB Graphics 

SCREEN 0 

CLS 

LOCATE 1, 1 
COLOR 9, 0 

PRINT "Do You Wish To Enter Data For Another Prediction (y/n)? 

COLOR 12, 0 
INPUT Answer? 

COLOR 15, 0 

Answer? = LCASE? (Answer? ) 

IF Answer? = "y" THEN GOTO PredictValue 
END 

f 

' End of main program. 

r 

GetMaterialProp: 'This subroutine searches for material property data records 

'in the material data file. 

GetMaterialPropl : 

OPEN MaterialDataFile? FOR INPUT AS #2 
FOR mat% = 1 TO NumMatDat% 

INPUT #2, dummy? 

NEXT mat% 

DO WHILE NOT EOF(2) 

INPUT #2, dummy? 

INPUT #2, TestMateriallD? 

IF TestMateriallD? = MateriallD? THEN 
FOR mat% = 1 TO NumMatDat% 

INPUT #2, Material# (mat%) 
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NEXT mat% 

CLOSE #2 
RETURN 
ELSE 

FOR mat% - 1 TO NumMatDat% + 1 
INPUT #2, dummy $ 

NEXT mat% 

END IF 
LOOP 

IF row% = -1 THEN 
SCREEN 0 
VIEW PRINT 
COLOR 15, 0 
CLS 

LOCATE 3, 1 

PRINT "ERROR - File ", MLITestDataFile$ 

PRINT "References Material Name " , MaterialID$ 

PRINT "That Is Not Contained In File ", MaterialDataFile$ 

PRINT 

PRINT "press any key to stop"; 

DO 

LOOP WHILE INKEY$ » " " 'Press any key to continue 

END 
END IF 
CLOSE #2 
COLOR 12, 0 
LOCATE row% + 1, 1 

PRINT "Material Name Not Found Please Re-enter (or enter QUIT to stop) 
LOCATE row%, col% 

PRINT " 

LOCATE row%, col% 

INPUT Material ID$ 

MaterialID$ = UCASE$(MaterialID$) 

IF Material ID$ = "QUIT" THEN END 
LOCATE row% + 1, 1 
PRINT " 

GOTO GetMaterialPropl 
RETURN 

f 

' This subroutine graphically displays the results. 

/ 

Graphics : 

SCREEN 9 
COLOR 15, 1 
CLS 

' Find largest function value and smallest and largest x values 
Fmaxl = PREDICTION# 

Xmin 1 « PredictVectMag# 

Xmax! ■ PredictVectMag# 

FOR it = 1 TO 10 

IF Fcalc# ( i% ) > Fmaxl THEN Fmaxl = Fcalc#(i%) 

IF Xcalc#(i%) < Xmin I THEN Xmin l = Xcalc#(i%) 

IF Xcalc#(i%) > Xmax l THEN Xmax l = Xcalc#(i%) 

NEXT i% 

Dxl = Xmax l - Xmin l 
wl = Dxl * .01 

WINDOW (-Dxl * .7, -Fmaxl * .17) -(Dxl * 1.2, Fmaxl * 1.2) 

LINE (-wl * 5, 0 ) - ( Dx 1 * 1.1, 0), 7 

LINE (-wl * 5, 0 ) - ( -w l * 5, Fmaxl * 1.1), 7 

LINE (-Dxl * .If -Fmaxl * .07) -(Dxl * 1.16, Fmaxl * 1.16), 7, B 
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FOR i% = 1 TO 10 

LINE (Xcalc#(i%) - Xminl - wl, 0)-(Xcalc#(i%) - Xminl + wt, Fcalc#(i%)), 11, 

NEXT i% 

/ 

'Display Interpolated Function 

COLOR 14, 1 

FOR i% = 1 TO 1001 

x# - Xminl + (i% - 1) * Dx! / 1000# 

f# = B#(l) + B# ( 2 ) * x# + B#(3) * x# * 2 + B# ( 4 ) * x# * 3 + B#(5) * x# * 4 
IF f# < 0 THEN f# = 0 
PSET (x# - Xminl, f#) 

NEXT i% 


LINE ( PredictVectMag# - Xminl - wl , -Fmaxl * . 005)-(PredictVectMag# - Xminl + wl 
LOCATE 2, 30 
COLOR 7, 1 
PRINT "FUNCTION" 

LOCATE 3, 30 
PRINT "MAGNITUDE" 

LOCATE 23, 56 

PRINT ” DISTANCE FROM ORIGIN 
COLOR 11, 1 
LOCATE 3, 2 

PRINT " BARS INDICATE" 

LOCATE 4, 2 

PRINT "DATABASE INTERPOLATIONS" 

LINE ( — Dxl * .68, Fmaxl * 1.05)-(-Dxl * .58, Fmaxl * 1.08), 11, BF 
LOCATE 7, 2 
COLOR 12, 1 

PRINT " BAR INDICATES" 

LOCATE 8, 2 

PRINT "FUNCTION PREDICTION" 

LINE (-Dxl * .68, Fmaxl * .83)-(-Dxl * .58, Fmaxl * .86), 12, BF 
LOCATE 11, 2 
COLOR 14, 1 

PRINT ” LINE INDICATES" 

LOCATE 12, 2 

PRINT "LEAST SQUARES FIT" 

LOCATE 13, 2 

PRINT "THROUGH INTERPOLATIONS" 

LINE (-Dxl * .68, Fmaxl * .625)-(-Dxl * .58, Fmaxl * .63), 14, BF 
LOCATE 20, 2 
COLOR 7, 1 

PRINT "PREDICTED VALUE:" 

LOCATE 21, 1 

IF Predict ionType% = 1 THEN 
PRINT " Bumper Hole" 

PRINT " Major Axis (in) = 

END IF 

IF PredictionType% = 2 THEN 
PRINT " Bumper Hole" 

PRINT " Minor Axis (in) = 

END IF 

IF PredictionType% = 3 THEN 
PRINT " MLI Hole" 

PRINT " Diameter (in) = 

END IF 

IF PredictionType% = 4 THEN 
PRINT " MLI Mass" 
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PRINT " Loss (grams) =" 

END IF 

IF PredictionType% • 5 THEN 
PRINT " Pressure Wall Hole" 

PRINT " Major Axis (in) =" 

END IF 

IF PredictionType% = 6 THEN 
PRINT " Pressure Wall Hole" 

PRINT " Minor Axis (in) = " 

END IF 

LOCATE 23, 4 
COLOR 15, 1 

PRINT USING "##.### PREDICTION# 

LOCATE 25, 28 
COLOR 12, 1 

PRINT "press any key to continue"; 

DO 

LOOP WHILE INKEY$ = "" 'Press any key to continue 

RETURN 

f 

TestDataFileError : 

COLOR 12, 0 
LOCATE 4, 1 

PRINT "Please Re-enter File Name (or enter QUIT to stop)" 
LOCATE row%, col% 

PRINT " 

LOCATE row%, col% 

INPUT MLITestDataFile$ 

MLITestDataFile$ - UCASE$ (MLITestDataFile$ ) 

IF MLITestDataFile$ = "QUIT" THEN END 
LOCATE 4, 1 
PRINT " 

RESUME 

9 

MaterialDataFileError: 

COLOR 12, 0 
LOCATE 4, 1 

PRINT "Please Re-enter File Name (or enter QUIT) to stop" 
LOCATE row%, col% 

PRINT " 

LOCATE row%, col% 

INPUT MaterialDataFile$ 

MaterialDataFile$ = UCASE$(MaterialDataFile$) 

IF MaterialDataFile$ = "QUIT" THEN END 
LOCATE 4, 1 
PRINT " 

RESUME 


( 
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POLYMETH . BAS 

Source code for the polynomial functions damage prediction program 


ECLARE SUB CalcFullAraatrix (i%, a#(), dx#()) 

ECLARE SUB CalcSmallAmatrix (i%, a#(), dx#()) 

ECLARE SUB CalcSmallestAmatrix (i%, a#(), dx#()) 

This program fits a linear polynomial through sets of 16 data records. 

The polynomial is of the form: 

DAMAGE - dl + d2*xl + d3*x2 + d4*x3 + d5*x4 

+ d6*xl*x2 + d7*xl*x3 + d8*xl*x4 + d9*x2*x3 + dl0*x2*x4 + dll*x3*x4 
+ dl2*xl*x2*x3 + dl3*xl*x2*x4 + dl4*xl*x3*x4 + dl5*x2*x3*x4 
+ dl6*xl*x2*x3*x4 

where: di = i-th coefficient 
xl = bumper thickness 
x2 = projectile diameter 
x3 = impact angle 
x4 = projectile velocity 

If not enough data records are available to fit the full function then 
the program (either automatically or as directed by the user by pressing 
function keys) will attempt to fit two lower order (also 
incomplete) polynomials of the following form: 

"simple": Damage = dl + d2*xl + d3*x2 + d4*x3 + d5*x4 

+ d6*xl*x2 + d7*xl*x3 + d8*xl*x4 + d9*x2*x3 + dl0*x2*x4 

"simplest" : Daamage = dl + d2*xl + d3*x2 + d4*x3 + d5*x4 
This program assumes that xl to x4 are the only parameters that vary signific 

Here the DAMAGE consists of bumper hole major and minor axis, 

MLI hole diameter and mass loss, and the pressure wall hole 
major and minor axis. 

lECLARE SUB weightedprediction ( ) 

> ECLARE SUB ParamWarning ( ) 

Dimensioned for 100 data points. 

:OMMON SHARED prediction# () , meandi stance# () , numsolutions% , f inalprediction# , n 
)IM activedof%(l TO 4) 'This variable keeps track of which of xl to x4 vary in t 
>IM TestID$(l TO 100) 'Data record ID. 

)IM DataSource$ ( 1 TO 100) 'Data record source. 

)IM TestDate$(l TO 100) 'Date of data record test. 

)IM BumperThickness#(l TO 100) 

)IM Pro jectileDiameter# ( 1 TO 100) 

)IM ImpactAngle# ( 1 TO 100) 

)IM Pro jectileVelocity# ( 1 TO 100) 

)IM BumperMajoraxis#( 1 TO 100) 

)IM BumperMinorAxis#(l TO 100) 

JIM MLIHoleDiam# ( 1 TO 100) 

JIM MLIMassLoss# ( 1 TO 100) 

JIM PressWallMa jAxis# ( 1 TO 100) 

JIM PressWallMinAxis#( 1 TO 100) 

JIM meandistance# ( 1 TO 10) 'Mean distance from prediction point to data points 
JIM datasort%(l TO 100) 'Vector to use for data sorting. 

JIM dx#(l TO 4) 'Distance from xi of data record to xi of prediction point, wh 
JIM a#(l TO 16, 1 TO 16) 'Matrix for solving for function coefficients: [a]{b} 
JIM b#(l TO 16) 'Function coefficients. 
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DIM c# ( 1 TO 16) 
DIM prediction# ( 1 
seed% = ((TIMER * 
RANDOMIZE seed% 
numsolutions% = 5 
CALL ParamWarning 


'Measured damage at a data point. 

«P0 5) 'The predictions from 5 sets of coefficients are stoi 
65536) / 86400) - 32768 'Used for generating random numbers. 

'Number of function fits to be calculated. 

'This subroutine warns the user that only four variables are \ 


CLS 

COLOR 9, 0 
LOCATE 1, 1 
PRINT CHR$( 201 ) ; 

FOR i% = 1 TO 78 
PRINT CHR$ ( 205 ) ; 

NEXT i% 

PRINT CHR$ (187) 

LOCATE 2, 1: PRINT CHR$(186) 

LOCATE 2, 80: PRINT CHR$(186) 

LOCATE 3, 1 
PRINT CHR$ ( 200 ) ; 

FOR i% - 1 TO 78 
PRINT CHR$ ( 205 ) ; 

NEXT i% 

PRINT CHR$ (188) , . . „ 

ON ERROR GOTO TestDataFileError 'This traps database file name problems. 

LOCATE 2, 3 
row% = 2 
col% - 24 
COLOR 11, 0 

PRINT "Test Data File Name? 

COLOR 12, 0 

INPUT MLITestDataFile$ 

OPEN MLITestDataFile$ FOR INPUT AS #1 
ON ERROR GOTO 0 

numdata% = 0 'Number of data records in the database. 

' The following variables store the average values of database data. 

BumpThkAve# = 0 

ProjDiaAve# = 0 

ProjVelAve# = 0 

ImpAngAve# = 0 

BumpMa jAxisAve# = 0 

BumpMinAxisAve# “ 0 

MLIHoleDiamAve# * 0 

MLIMassLossAve# * 0 

PressWallMa jAxisAve# =» 0 

PressWallMinAxisAve# = 0 

'Modify top line of boarder. 

COLOR 9, 0 
LOCATE 3, 1 
PRINT CHR$ (204 ) 

LOCATE 3, 80 
PRINT CHR$ ( 185 ) 

VIEW PRINT 4 TO 12 

' Read in data from database file. 


DO WHILE NOT EOF(l) 

numdata% = numdata% + 1 
INPUT #1, dummy $ 

INPUT #1, TestID$(numdata%) 

INPUT #1, DataSource$(numdata%) 

INPUT #1, TestDate$(numdata%) 

COLOR 9, 0 160 

PRINT CHR$ ( 186 ) ; 
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COLOR 11, 0 

' Scroll data to screen as it is read. 

PRINT " No . : " ; 

PRINT numdata%; 

PRINT " ID: " ; 

PRINT Test ID$ ( numdata% ) ; 

PRINT " Source: 

PRINT DataSource$ (numdata% ) ; 

PRINT " Date: 

PRINT TestDate$ ( numdata% ) ; 

LOCATE CSRLIN, 80 
COLOR 9, 0 
PRINT CHR$ ( 186 ) 

' Input bumper material properties . 

INPUT #1, dummy $ 

INPUT #1, BumperThickness#(numdata%) 

BumpThkAve# = BumpThkAve# + BumperThickness#(numdata%) 

INPUT #1, dummy# 

INPUT #1, dummy $ 

INPUT #1, dummy# 

INPUT #1, dummy $ 

INPUT #1, ProjectileDiameter#(numdata%) 

ProjDiaAve# = ProjDiaAve# + Pro jectileDiameter#(numdata% ) 

INPUT #1, ImpactAngle#(numdata%) 

ImpAngAve# = ImpAngAve# + ImpactAngle#(numdata%) 

INPUT #1, ProjectileVelocity#(numdata%) 

ProjVelAve# = ProjVelAve# + Pro jectileVelocity#(numdata% ) 

INPUT #1, BumperMajoraxis#(numdata%) 

BumpMa jAxisAve# ■ BumpMa jAxisAve# + BumperMa joraxis#(numdata% ) 

INPUT #1, BumperMinorAxis# ( numdata% ) 

BumpMinAxisAve# = BumpMinAxisAve# + BumperMinorAxis# (numdata% ) 

INPUT #1, MLIHoleDiam#(numdata%) 

MLIHoleDiamAve# = MLIHoleDiamAve# + MLIHoleDiam#(numdata%) 

INPUT #1, MLIMassLoss# (numdata% ) 

MLIMassLossAve# = MLIMassLossAve# + MLIMassLoss# (numdata% ) 

INPUT #1, PressWallMajAxis#(numdata%) 

PressWallMa jAxisAve# = PressWallMa jAxisAve# + PressWallMa jAxis# (numdata% ) 
INPUT #1, PressWallMinAxis# ( numdata% ) 

PressWallMinAxisAve# =* PressWallMinAxisAve# + PressWallMinAxis# (numdata% ) 
INPUT #1, dummy$ 

LOOP 

VIEW PRINT 


'build box for averages 
LOCATE 12, 2 
FOR i% * 1 TO 78 
PRINT CHR$ ( 205 ) ; 

NEXT i% 

LOCATE 12, 1: PRINT CHR$(204) 

LOCATE 12, 40: PRINT CHR$(203) 

LOCATE 12, 80: PRINT CHR$(185) 

FOR i% = 13 TO 23 

LOCATE i% , 1: PRINT CHR$(186) 

LOCATE i% , 40: PRINT CHR$(186) 

LOCATE i%, 80: PRINT CHR$(186) 

NEXT i% 

LOCATE 17, 1: PRINT CHR$(204); 

FOR i% = 1 TO 78 

IF i% <> 39 THEN PRINT CHR$(205); 161 

IF i% = 39 THEN PRINT CHR$(206); 



POLYMETH.BAS 


NEXT i% 

PRINT CHR$(185) 

LOCATE 24, 1: PRINT CHR$(200); 

FOR i% = 1 TO 78 

IF i% <> 39 THEN PRINT CHR$(205); 

IF i% = 39 THEN PRINT CHR$(202); 

NEXT i% 

' Print data averages to the screen. 

PRINT CHR$ (188) 

BumpThkAve# = BumpThkAve# / numdata% 

LOCATE 13, 3 
COLOR 11, 0 

PRINT "Ave . Bumper Thk ( in ) : " ; 

COLOR 15, 0 

PRINT USING "##.### A AAA " ; BumpThkAve#; 

ProjDiaAve# = ProjDiaAve# / numdata% 

LOCATE 13, 42 
COLOR 11, 0 

PRINT "Ave. Proj. Dia. (in):"; 

COLOR 15, 0 

PRINT USING "##.### "; ProjDiaAve#; 

ImpAngAve# = ImpAngAve# / numdata% 

LOCATE 15, 3 
COLOR 11, 0 

PRINT "Ave* Impact Angle (deg):"; 

COLOR 15, 0 

PRINT USING "##.### "; ImpAngAve# ; 

ProjVelAve# = ProjVelAve# / numdata% 

LOCATE 15, 42 
COLOR 11, 0 

PRINT "Ave. Proj. Vel. (km/sec):"; 

COLOR 15, 0 

PRINT USING "##.### "; ProjVelAve#; 

BumpMa jAxisAve# = BumpMa jAxisAve# / numdata% 

LOCATE 19, 3 
COLOR 11, 0 

PRINT "Ave. Maj. Bumper Hole (in):"; 

COLOR 15, 0 

PRINT USING "##.### AAAA " ; BumpMa jAxisAve#; 
BumpMinAxisAve# = BumpMinAxisAve# / numdata% 

LOCATE 19, 42 
COLOR 11, 0 

PRINT "Ave. Min. Bumper Hole (in):"; 

COLOR 15, 0 

PRINT USING "##.### A "; BumpMinAxisAve# ; 

MLIHoleDiamAve# = MLIHoleDiamAve# / numdata% 

LOCATE 21, 3 
COLOR 11, 0 

PRINT "Ave. MLI Hole Diam. (in):"; 

COLOR 15, 0 

PRINT USING "##.### AAAA " ; MLIHoleDiamAve#; 
MLIMassLossAve# = MLIMass Loss Ave# / numdata% 

LOCATE 21, 42 
COLOR 11, 0 

PRINT "Ave. MLI Mass Loss (grams):"; 

COLOR 15, 0 

PRINT USING "##.### AAAA " ; MLIMassLossAve#; 

PresBWallMa jAxisAve# = Pres sWallMa jAxisAve# / numdata% 
LOCATE 23, 3 
COLOR 11, 0 
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UNT "Ave. Maj. P.Wall Hole (in)s"; 

DLOR 15, 0 

UNT USING "##.### ~ A A ; PressWallMa jAxisAve# ; 
ressWallMinAxisAve# ■ PressWallMinAxisAve# / numdata% 
DCATE 23, 42 
DLOR 11, 0 

UNT "Ave. Min. P.Wall Hole (in)j H ; 

DLOR 15, 0 

UNT USING ”##.### PressWallMinAxisAve#? 

DCATE 25, 28 
DLOR 12, 0 

UNT "press any key to continue"; 

DEEP 


Check which DOF are "active" - ie. vary in the database. 

If DOF xi is active then activedof%(i) is set equal to 1 (equal to 0 if not a 
Also, normalize the data wrt the average. 

DR i% = 1 TO numdata% 

BumperThickness# ( i% ) = BumperThickness# (i% ) / BumpThkAve# 

IF BumperThickness# ( i% ) < .999 OR BumperThickness# (i%) > 1.001 THEN activedof 
Pro jectileDiameter# ( i% ) = ProjectileDiameter#(i%) / ProjDiaAve# 

IF Pro jectileDiameter# ( i% ) < .999 OR Pro jectileDiameter# ( i% ) > 1.001 THEN act 
IF ImpAngAve# <> 0 THEN 

ImpactAngle# ( i% ) = ImpactAngle# ( i% ) / ImpAngAve# 

IF ImpactAngle# (i%) < .999 OR ImpactAngle# (i%) > 1.001 THEN activedof %( 3 ) 
END IF 

Pro jectileVelocity# ( i% ) = ProjectileVelocity#(i%) / ProjVelAve# 

IF Pro jectileVelocity# (i% ) < .999 OR Pro jectileVelocity# ( i% ) > 1.001 THEN act 
2XT i% 

imactivedof % = 0 'This is the total number of active DOF (can range from 1 to 4 
DR i% = 1 TO 4 

numactivedof % = numactivedof % + activedof % (i%) 

2XT i% 


illFuncReqRecords% = 2 ~ numactivedof % 'Calculate number of records required 

? numactivedof % = 1 THEN ReqRecords% = 2 'Calculate number of records requir 

? numactivedof % = 2 THEN ReqRecords% = 4 

? numactivedof % = 3 THEN ReqRecords% = 7 

? numactivedof % = 4 THEN ReqRecords% = 11 

LnReqRecords% = 1 + numactivedof % 'Find smallest allowable number of records r 
f FullFuncReqRecords% > numdata% THEN 'Issue warning and stop if not enough dat 


CLS 

COLOR 12, 0 

PRINT "Warning - at least " ; 

COLOR 11, 0 

PRINT FullFuncReqRecords%; 

COLOR 12, 0 

PRINT " data records are required to fit the full function!" 
PRINT 

PRINT "Data file "; 

COLOR 11, 0 

PRINT MLITestDataFile$ ? 

COLOR 12, 0 
PRINT " only has " ; 

COLOR 11, 0 
PRINT numdata%; 

COLOR 12, 0 
PRINT " records . " 

IF numdata% >= MinReqRecords% THEN 
COLOR 13, 0 
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PRINT "Fitting incomplete polynomial will now be attempted l" 

PRINT 

PRINT 

COLOR 14, 0 

PRINT "Press any key to continue." 

SLEEP 

GOTO SetDefaultValues 
END IF 
PRINT 
PRINT 

COLOR 14, 0 

PRINT "Press any key to stop." 

SLEEP 
STOP 
END IF 

SetDefaultValues: 'Set the default values - these can be changed. 
PredictBumpThick# = .04 / BumpThkAve# 

PredictPro jDia# ■ .313 / ProjDiaAve# 

IF ImpAngAve# <> 0 THEN Predict ImpAngle# = 45 / ImpAngAve# 

PredictPro jVel# = 5.3# / ProjVelAve# 

' Top of prediction loop, 
predictvalue : 

t 

' Here the default order of the function is determined. Variable amatrixcode% k 
IF numdata% >= NumRequiredDataRecord% THEN amatrixcode% = 0 'Attempt to fit ful 
IF numdata% >= ReqRecords% AND numdata% < NumRequiredDataRecord% THEN amatrixcod 
IF numdata% >= MinReqRecords% AND numdata% < ReqRecords% THEN amatrixcode% = 2 ' 

t 

COLOR 15, 1 
CLS 

' Prompt user for the required prediction. 

LOCATE 1, 3 

PRINT "Please Enter The ■ ; 

COLOR 12, 1 

PRINT "Number " ; ' 

COLOR 15, 1 

PRINT "Associated With the Desired Action:" 

LOCATE 5, 3 

COLOR 12, 1 

PRINT "1. "; 

COLOR 15, 1 

PRINT "Predict Bumper Hole Major Axis" 

LOCATE 7, 3 ~ 

COLOR 12, 1 

PRINT "2. "; -- 

COLOR 15, 1 ' 

PRINT "Predict Bumper Hole Minor Axis" 

LOCATE 9, 3 

COLOR 12, 1 - 

PRINT "3. "? 

COLOR 15, 1 

PRINT "Predict MLI Hole Diameter" 

LOCATE 11, 3 
COLOR 12, 1 

PRINT "4. "? 

COLOR 15, 1 

PRINT "Predict MLI Mass Loss" 

LOCATE 13, 3 i 164 

COLOR 12, 1 

PRINT "5. "; 
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OLOR 15 , 1 

RINT "Predict Pressure Wall Hole Major Axis" 

OCATE 15, 3 
OLOR 12, 1 
RINT "6. 

OLOR 15, 1 

RINT "Predict Pressure Wall Hole Minor Axis" 

OCATE 17, 3 
OLOR 12, 1 
RINT "7. 

OLOR 15, 1 
RINT "Quit" 
electPredictionType : 

.OCATE 19, 1 
lOLOR 12, 1 

INPUT "Enter Number" ; predictiontype% 

;F predictiontype% < 1 OR predictiontype% > 7 THEN 
LOCATE 19, 1 
PRINT ” 

GOTO SelectPredictionType 
1ND IF 

:f predictiontype% = 7 THEN END 
'.OLOR 11, 0 

:ls 

.OCATE 1, 1 

’RINT "ENTER DATA FOR DESIRED PREDICTION: " 

IOLOR 15, 0 
.OCATE 3, 1 

’RINT "[default values shown in square brackets]" 

IOLOR 10, 0 
.OCATE 5, 1 

’RINT "(magnitude relative to database average shown in round brackets)" 
)OFWaming% = 0 
.OCATE 8, 1 

’redictBumpThick# = PredictBumpThick# * BumpThkAve# 

IOLOR 11, 0 

’RINT "Bumper Thickness (in):"; 
iolor 15, 0 
’RINT "["; 

’RINT USING "## • ###*^ /v " ; PredictBumpThick#; 

’RINT "]" 

.OCATE 8, 40 
IOLOR 12, 0 
INPUT "", dummy# 

IF dummy# <> 0 THEN PredictBumpThick# = dummy# 

’redictBumpThick# ■ PredictBumpThick# / BumpThkAve# 

IF (PredictBumpThick# < .999 OR PredictBumpThick# > 1.001) AND activedof % ( 1 ) 
.OCATE 8, 60 
IOLOR 10, 0 
’RINT "("; 

’RINT USING "##.###" ; PredictBumpThick#; 

’RINT " )" 

[F DOFWarning% = 1 THEN 
LOCATE 10, 1 
COLOR 12, 0 

PRINT "Waming-Your test data does not support the above input data 1 " 

2ND IF 

)OFWaming% = 0 
I.OCATE 12, 1 

’redictProjDia# = PredictPro jDia# * ProjDiaAve# 
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COLOR 11, 0 

PRINT "Projectile Diameter (in)*"? 

COLOR 15, 0 
PRINT "["? 

PRINT USING "##.### AAAA "? PredictPro jDia# ? 

PRINT "]" 

LOCATE 12, 40 
COLOR 12, 0 
INPUT "", dummy# 

IF dummy# <> 0 THEN PredictPro jDia# = dummy# 

PredictPro jDia# ■ PredictPro jDia# / ProjDiaAve# 

IF (PredictPro jDia# < .999 OR PredictPro jDia# > 1.001) AND activedof%(2) 
LOCATE 12, 60 
COLOR 10, 0 
PRINT "("? 

PRINT USING "##.###"? PredictPro jDia#? 

PRINT ")" 

IF DOFWarning% = 1 THEN 


LOCATE 14, 1 
COLOR 12, 0 
PRINT "Warning 
END IF 


-Your test data does not support the above input data l " 


0 THE 


DOFWaming% = 0 
LOCATE 16, 1 
COLOR 11, 0 

PRINT "Impact Angle (degrees):"? 

LOCATE 16, 40 
COLOR 12, 0 

INPUT "", Predict ImpAngle# 

IF ImpAngAve# <> 0 THEN 

PredictlmpAngle# ■ Predict ImpAngle# / ImpAngAve# 

IF (PredictlmpAngle# < .999 OR PredictlmpAngle# > 1.001) AND activedof%(3) - 
LOCATE 16, 60 
COLOR 10, 0 
PRINT "("? 

PRINT USING "##.###"? PredictlmpAngle#? 

PRINT ")" 

IF DOFWaming% = 1 THEN 


LOCATE 18, 1 
COLOR 12, 0 

PRINT "Warning-Your test data does not support the above input datal 
END IF 
END IF 

IF ImpAngAve# = 0 AND PredictlmpAngle# > 0 THEN 


LOCATE 18, 1 
COLOR 12, 0 

PRINT "Warning- Your test data does not support the above input data J " 
END IF 


H 


DOFWarning% = 0 
LOCATE 20, 1 

PredictPro jVel# = PredictPro jVel# * ProjVelAve# 
COLOR 11, 0 

PRINT "Proj. Vel . (km/sec) :"? 

COLOR 15, 0 
PRINT "["? 

PRINT USING "##.### AAAA " ? PredictPro jVel#? 

PRINT "]"? 

LOCATE 20, 40 166 

COLOR 12, 0 

INPUT ? dummy# 
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IF dummy# <> 0 THEN PredictPro jVel# => dummy# 

PredictPro jVel# = PredictPro jVel# / ProjVelAve# 

IF (PredictProjVel# < .999 OR PredictPro jVel# > 1.001) AND activedof % ( 4 ) = 0 THE 
LOCATE 20, 60 
COLOR 10, 0 
PRINT 

PRINT USING "##.###"; PredictProjVel#; 

PRINT " 

IF DOFWarning% = 1 THEN 
LOCATE 22, 1 
COLOR 12, 0 

PRINT "Warning-Your test data does not support the above input data ! " ; 

END IF 

LOCATE 25, 28 
COLOR 12, 0 

PRINT "press any key to continue"; 

SLEEP 

CLS 

calculationattempts% = 0 'Stores the number of attempts made to seek function 
IF amatrixcode% = 0 THEN 

LOCATE 21, 20: COLOR 10, 0 

PRINT "Currently Seeking Full Function" 

END IF 

IF amatrixcode% = 1 THEN 

LOCATE 21, 20: COLOR 10, 0 

PRINT "Currently Seeking Simpler Function" 

END IF 

IF amatrixcode% = 2 THEN 

LOCATE 21, 20: COLOR 10, 0 

PRINT "Currently Seeking Simplest Function" 

END IF 

LOCATE 22, 22: COLOR 12, 1 

PRINT " PRESS FI TO STOP CALCULATIONS 

LOCATE 23, 22: COLOR 15, 1 

PRINT " PRESS F2 TO TRY SIMPLER FUNCTION "; 

LOCATE 24, 22: COLOR 11, 1 

PRINT " PRESS F3 TO TRY SIMPLEST FUNCTION "; 
solution% = 1 'Stores the number of solutions obtained. 

FindSolution : 

' Use function keys to allow user to select the function to be fit. 

ON KEY(l) GOSUB stopcalculations 'User can stop calculations if no solutions 
KEY ( 1 ) ON 

ON KEY(2) GOSUB simplifycalculations 'User can request simpler function if no 
KEY ( 2 ) ON 

ON KEY(3) GOSUB simplestcalculations 'User can request simpler function if no 
KEY ( 3 ) ON 

InitialDataSort : 'Here we randomly select data points for the function fit. 

FOR i% = 1 TO 2 * numactivedof % 
tryagainl : 

try% = numdata% * RND + 1 

IF try% > numdata% THEN try% = numdata% 

IF i% > 1 THEN 

FOR j% = 1 TO i% - 1 

IF try% = datasort% ( j% ) THEN GOTO tryagainl 'Make sure the same data 
NEXT j% 

datasort%(i% ) = try% 

ELSE 

datasort% ( 1 ) = try% 

END IF 
NEXT i% 
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0 # 

OTHEN order% = FullFuncReqRecords% 'Full function. 

1 THEN order% = ReqRecords% 'Simpler function. 

- 'Simplest function. 


ON ERROR GOTO RankError 'If selected data points are not linearly independer 

LOCATE 20, 20$ COLOR 14, 0 

PRINT "Working on calculation attempts "? 

LOCATE 20, 55s COLOR 15, 0 
PRINT calculationattempts% 

CalculateCoef f icients s 

calculationa ttempt s % - calculationattempts% + 1 

LOCATE 20, 55s COLOR 15, 0 
PRINT calculationattempts% 

COLOR 15, 0 

meandistance# ( solution% ) 

IF amatrixcode% = 

IF amatrixcode% ■ * . - „ , „ 

IF amatrixcode% = 2 THEN order% = MinReqRecords% 

FOR i% = 1 TO order % 

FOR j% = 1 TO 4 
dx#(j%) = 0# 

j% X = l % 'This counter is used to ensure that only active DOF are placed : 

dx#(j%) = BumperThickness#(datasort%(i%) ) - PredictBumpThick# 

j% = j% + 1 
END IF 

IF dJ#tj%t 0 ” % ProjectileDiameter#(datasort%(i%)) - PredictPro jDia# 

j% = j% + 1 
END IF 

IF activedof % ( 3 ) = 1 THEN . - 

dx#(j%) = ImpactAngle#(datasort%(i%) ) - PredictlmpAngle# 

j% - j% + 1 
END IF 

IF activedof %( 4 ) ■ 1 THEN , 

dx#(j%) - Pro jectileVelocity# (datasort% ( i% ) ) ~ PredictPro jVe # 

END IF 

distance# m 0# 

FOR j% = 1 TO 4 

distance# ■ distance# + dx#(j%) A 2 

distance# = SQR ( distance# ) 'Calculate the "distance" in the design space 
meandistance# (solution%) = meandistance# (solution%)+ distance# 

IF amatrixcode% = 0 THEN CALL CalcFullAmatrix(i%, a#(), dx#( ) 

IF amatrixcode% = 1 THEN CALL CalcSmallAmatrix(i%, a#(), dx#()) 


2 THEN CALL CalcSmallestAmatrix( i% , a#( ) f dxf ( ) ) 

• 3 - ft / 1 ft 


IF amatrixcode% = 

IF predictiontype% 

IF predictiontype% - 
IF predictiontype% * 

IF predictiontype% = 

IF predictiontype% * 

IF predict iontype% = 

NEXT i% 

meandistance# ( solution% ) 

The following statements are a Gauss Elimination Solver. 


THEN c#(i%) 
THEN c#(i%) 
THEN c# ( i% ) 
THEN c#(i%) 
THEN c#(i%) 
THEN c#(i%) 


liico uruiio y j. w / \ ; r \ t t 

BumperMa joraxis# (datasort% ( i% ) ) 
BumperMinorAxis# ( datasort% ( i% ) ) 
MLIHoleDiam# ( datasort % ( i% ) ) 
MLIMassLoss# (datasort% ( i% ) ) 
PressWallMa jAxis# (datasort% ( i% ) ) 

PressWallMinAxis# (datasort% ( i% ) ) 


meandistance# (solution%) / order% 


FOR i% = 1 TO order% - 1 
FOR j% = i% + 1 TO order% 

IF a#(j%, i%) <= a#(i%, i%) THEN GOTO LI 
FOR k% = i% TO order% 

dumb# = a#(i%, k%) 168 

a# (i% , k%) = a#(j%, k%) 
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a#(j%, k%) = dumb# 

NEXT k% 
dumb# = c#(i%) 
c#(i%) = c# ( j % ) 
c# ( j % ) = dumb# 

Lis 

NEXT j% 

FOR j% = i% + 1 TO order% 
factor# = a# ( j% , i%) / a#(i%, i%) 

FOR k% = i% TO order% 

a#(j%, k%) = a#(j%, k% ) - a#(i!, k% ) * factor# 

NEXT k% 

c#(j%) = c#(j%) - c#(i%) * factor# 

NEXT j! 

NEXT i% 

FOR i% = order % TO 1 STEP -1 
sum# = 0 

IF i% = order % THEN 

b#(i%) = (c#(i%) - sum#) / a#(i%, i%) 

GOTO L2 
END IF 

FOR j! = i! + 1 TO order% 

sum# = sum# + a#(i%, j%) * b#(j%) 

NEXT j% 

b#(i%) * (c#(i%) - sum#) / a#(i!, i%) 

L2 : 

NEXT i% 

9 

' End of Gauss Elimination Solver. 

9 

prediction# ( solution% ) = b#(l) 'The constant coefficient is the prediction - 
' Here a check is made to ensure that the calculated result is reasonable. 

IF predictiontype% = 1 THEN 

IF prediction# ( solution% ) < 0# OR prediction# ( solution% ) > 100# * BumpMajA 
END IF 

IF predictiontype% = 2 THEN 

IF prediction# ( solution% ) < 0# OR prediction# ( solution% ) > 100# * BumpMinA 
END IF 

IF predictiontype! = 3 THEN 

IF prediction# ( solution! ) < 0# OR prediction# ( solution! ) > 100# * MLIHoleD 
END IF 

IF predictiontype! = 4 THEN 

IF prediction# ( solution! ) < 0# OR prediction# (solution!) > 100# * MLIMassL 
END IF 

IF predictiontype! = 5 THEN 

IF prediction# ( solution! ) < 0# OR prediction# ( solution! ) > 100# * PressWal 
END IF 

IF predictiontype! - 6 THEN 

IF prediction# (solution!) < 0# OR prediction# ( solution! ) > 100# * PressWal 
END IF 
KEY ( 1 ) OFF 
KEY ( 2 ) OFF 
KEY ( 3 ) OFF 
COLOR 11, 0 

LOCATE 2 * solution! - 1, 1 

IF predictiontype! = 1 THEN PRINT "Predicted Bumper Hole Major Axis (in) 

IF predictiontype! = 2 THEN PRINT "Predicted Bumper Hole Minor Axis (in) ="; 

IF predictiontype! * 3 THEN PRINT "Predicted MLI Hole Diameter (in) ="; 

IF predictiontype! = 4 THEN PRINT "Predicted MLI Mass Loss (grams) ="; 

IF predictiontype! = 5 THEN PRINT "Pred. Press Wall Hole Major Axis (in) ="; 
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IF predictiontype% - 6 THEN PRINT "Pred. Press Wall Hole Minor Axis (in) 
COLOR 15, 0 

PRINT USING prediction# ( solution* ) ; 

COLOR 14, 0 

PRINT " mean data dist. - 
COLOR 15, 0 

PRINT USING " ## . ###~~ ~ ~ meandistance# ( solution*) 

solution* * solution* + 1 

IF solution* <= numsolutions* THEN GOTO FindSolution 
CALL weightedprediction 
LOCATE 13, 1 
COLOR 10, 0 

PRINT "Weighted Sum of Above Predictions: 

COLOR 15, 0 

PRINT USING "##.###'■''""; f inalprediction# 

LOCATE 20, 20 
PRINT ■ 

LOCATE 21, 1 
PRINT " 

LOCATE 22, 20 
PRINT " 

LOCATE 23, 20 
PRINT " 

LOCATE 24, 20 

PRINT ■ " ' 

LOCATE 23, 1 
COLOR 9, 0 

PRINT "Do you wish to make more predictions (y/n)? 

COLOR 12, 0 
INPUT answer$ 

COLOR 15, 0 
IF answer$ =* ■ ” THEN 
LOCATE 23, 1 
COLOR 9, 0 

PRINT "Do you wish to make more predictions (y/n)? "; 

COLOR 12, 0 
INPUT "", answers 
COLOR 15, 0 
END IF 

answers - LCASE$( answerS) 

COLOR 15, 0 

IF answers ■ "y" THEN GOTO predictvalue 
END 'End of program. 

r 

TestDataFileError: 

COLOR 12, 0 
LOCATE 4, 1 

PRINT "Please Re-enter File Name (or enter QUIT to stop)" 

LOCATE row*, col* 

PRINT " 

LOCATE row*, col* 

INPUT " " , MLITestDataFileS 

MLITestDataFileS » UCASE$ (MLITestDataFileS ) 

IF MLITestDataFileS - "QUIT" THEN END 
LOCATE 4, 1 
PRINT " 

RESUME 


RankError: 

FOR i* = 1 TO 2 ‘ numactivedof % 
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tryagain s 

try% = numdata% * RND + 1 

IF try% > numdata% THEN try% ■ numdata% 

IF i% > 1 THEN 

FOR j% - 1 TO i% - 1 

IF try% = datasort%( j%) THEN GOTO tryagain 
NEXT j% 

datasort% ( i% ) = try% 

ELSE 

datasort%(l) = try% 

END IF 
NEXT i% 

RESUME CalculateCoef f icients 

t 

stopcalculations : 

END 

t 

simplif ycalculations : 
amatrixcode% = 1 
LOCATE 21 , 1 
PRINT " 

LOCATE 21, 20 
COLOR 10, 0 

PRINT "Currently Seeking Simpler Function" 

RESUME CalculateCoef f icients 
/ 

simplestcalculations : 
amatrixcode% = 2 
LOCATE 21, 1 
PRINT " 

LOCATE 21, 20 
COLOR 10, 0 

PRINT "Currently Seeking Simplest Function" 

RESUME CalculateCoef f icients 

SUB CalcFullAmatrix (i%, a#(), dx#()) 

' Here the coefficients for the complete function are calculated. 

a#(i%, 1) = 1# 

a#(i%, 2) = dx# ( 1 ) 

a#(i%, 3) = dx# ( 2 ) 

a#(i%, 4) = dx# ( 1 ) * dx# ( 2 ) 

a#(i%, 5) = dx# ( 3 ) 

a#(i%, 6) = dx# ( 1 ) * dx# ( 3 ) 

a# ( i% , 7) = dx# ( 2 ) * dx#(3) 

a#(i%, 8) = dx# ( 1 ) * dx# ( 2 ) * dx#(3) 

a#(i%, 9) = dx# ( 4 ) 

a#(i%, 10) = dx# ( 1 ) * dx# ( 4 ) 

a# (i%, 11) = dx# ( 2 ) * dx# ( 4 ) 

a#(i%, 12) = dx# ( 3 ) * dx#(4) 

a# (i%, 13) = dx# ( 1 ) * dx#( 2 ) * dx#(4) 

a# (i%, 14) = dx# ( 1 ) * dx#( 3 ) * dx#(4) 

a#(i%, 15) = dx# ( 2 ) * dx#(3) * dx#(4) 

a#(i%, 16) = dx# ( 1 ) * dx# ( 2 ) * dx#(3) * dx#(4) 

END SUB 

SUB CalcSmallAmatrix (i%, a#(), dx#()) 

' Here the coefficients for the simpler function are evaluated. 
a#(i%, 1) = 1# 

a#(i%, 2) = dx# ( 1 ) 171 

a# ( i% , 3) = dx# ( 2 ) 
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a#(i%, 4) = dx# ( 1 ) * dx# ( 2 ) 
a#(i%, 5) = dx# ( 3 ) 
a#(i%, 6) = dx# ( 1 ) * dx# ( 3 ) 
a# (i%, 7) = dx# ( 2 ) * dx#(3) 
a#(i%, 8) = dx# ( 4 ) 
a# (i%, 9) = dx# ( 1 ) * dx# ( 4 ) 
a#(i%, 10) = dx# ( 2 ) * dx#(4) 
a#(i%, 11) - dx# ( 3 ) * dx# ( 4 ) 

END SUB 

SUB CalcSmallestAmatrix (i%, a#(), dx#()) 

• Here the coefficients for the simpler function are evaluated. 
a#(i%, 1) = 1# 

a#(i%, 2) = dx# ( 1 ) 
a#(i%, 3) = dx# ( 2 ) 
a#(i%, 4) = dx# ( 3 ) 
a#(i%, 5) = dx# ( 4 ) 

END SUB 

SUB ParamWarning 
COLOR 12, 0 
CLS 

LOCATE 1 1 

PRINT "WARNING - This program assumes that only the following system" 
PRINT " parameters vary significantly in the database*" 

PRINT 

COLOR 14, 0 

PRINT " 1* Bumper Thickness" 

PRINT " 2. Projectile Diameter" 

PRINT " 3. Projectile Velocity" 

PRINT " 4. Impact Angle" 

PRINT 

COLOR 12, 0 , ^ 

PRINT "All other system parameters are assumed to be constant througnout^ 
PRINT "the entire database or are assumed to have no influence on the" 
PRINT "amount of impact damage sustained." 

PRINT 

COLOR 9, 0 

PRINT "Do you wish to continue (y/n)? "; 

COLOR 12, 0 
INPUT " " , answers 

IF answers = " " THEN answers * "y" 
answers = LCASE$( answer$) 

IF answers = "n" THEN END 
COLOR 15, 0 
END SUB 

SUB weightedprediction 

* Here a weighted prediction is made based on all calculated results . 

' Find the two predictions with largest meandistance# values 

' and ignor these values when making the weighted prediction. 

largemeandistl# * meandistance# ( 1 ) 

numlargemeandistl% = 1 

FOR i% - 2 TO numsolutions% 

IF meandistance# (i%) > largemeandistl# THEN 
largemeandistl# » meandistance# (i%) 
numlargemeandistl% = i% 

END IF 172 

NEXT i% 

IF numlargemeandistl% <> 1 THEN 
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largemeandist2# * meandistance#( 1 ) 
numlargemeandistl% = 1 
ELSE 

largemeandist2# = meandi stance# ( 2 ) 
nvunlargemeandistl% ■ 2 
END IF 

FOR i% = 1 TO numsolutions% 

IF i% <> numlargemeandistl% THEN 

IF meandistance#(i%) > largemeandist2# THEN 
largemeandist2# = meandistance# ( i% ) 
numlargemeandist2% = i% 

END IF 
END IF 
NEXT i% 
factor# = 0# 
finalprediction# = 0# 
n% = nuinactivedof % - 1 
IF n% < 1 THEN n% = 1 
FOR it = 1 TO numsolutions% 

IF i% <> numlargemeandistl% AND i% <> numlargemeandist2% THEN 

finalprediction# = finalprediction# + prediction# ( i% ) / meandistance# (i%) 
factor# = factor# + 1# / meandistance# (i% ) ~ n% 

END IF 
NEXT i% 

finalprediction# = finalprediction# / factor# 

END SUB 
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NONDIMEN. BAS . . .. . . „ 

Source code for the nondimensional functions damage prediction program 

This program was written by William K. Rule, University of Alabama, (205)3 
This program makes predictions for the following functions : 

1. Bumper hole minor diameter. 

2. Bumper hole major diameter. 

3. ML I hole diameter. 

4. Pressure wall hole diameter. 

This program uses functions of the form given in the report* ___ 4 . DV.«n 

Schonberg, W. P., Bean, A. J., and Darzi, K. , M Hypervelocity Impact Phy 

DECLARE SUB OptParameters ( Inputlteration#, InputAlpha#) 

DECLARE SUB ShowCoef f icients (RSquaredValues#( ) ) ... 

DECLARE SUB Ob jectiveFunction (calc%, code%, xtry#(), objective#, rsquare #) 
DECLARE SUB DisplayConvergence (col%, rsquared#) 

DECLARE SUB RuleOpt (Inputlteration#, InputAlpha#, RSquared Values# ( ) ) 

DECLARE SUB MakePrediction ( ) 

COMMON SHARED PredictBumpThick# , PredictBumpStandOff#, PredictPressWallThick#, P 
COMMON SHARED pi#, a#(), BumperSoundSpeed# , nurodata% „ 

COMMON SHARED BumpMinDiaPred# , BumpMaxDiaPred# , MLIDiaPred#, p 

COMMON SHARED BumperThickness# ( ) , Bumper StandOf f#() , P ^®®^®"f^T^ Ck JJr IMa ^£ os 
COMMON SHARED BumperMa jorAxis# ( ) , BumperMinorAxis#( ) , MLIHoleDiam#( ) , 

$ 

• Vector a#() stores the function coefficients. 

DIM a# (1 TO 23) 

* Dimensioned for 100 data points. 

DIM TestID$(l TO 100) 

DIM DataSource$ ( 1 TO 100) 

DIM TestDate$ ( 1 TO 100) 

DIM BumperThickness# ( 1 TO 100) 

DIM BumperStandOff #( 1 TO 100) 

DIM pressureWallThickness# ( 1 TO 100) 

DIM ProjectileDiameter#(l TO 100) 

DIM ImpactAngle# ( 1 TO 100) 

DIM ProjectileVelocity#(l TO 100) 

DIM BumperMa jorAxis# (1 TO 100) 

DIM BumperMinorAxis#(l TO 100) 

DIM MLIHoleDiam# ( 1 TO 100) 

DIM MLIMassLoss# ( 1 TO 100) 

DIM PressWallMajAxis#(l TO 100) 

DIM PressWallMinAxis# ( 1 TO 100) .* on f or the 

' Vector RSquaredValues# ( ) stores the coefficients of determination for the 

DIM RSquaredValues# ( 4 ) 

9 

pi# = 3.14159265359# 

r 

COLOR 9, 0 
CLS 

LOCATE 1, 1 
PRINT CHR$(201) ; 

FOR i% = 1 TO 78 
PRINT CHR$ ( 205 ) ; 

NEXT i% 174 

PRINT CHR$ ( 187 ) 

LOCATE 2, Is PRINT CHR$(186) 
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.OCATE 2, 80: PRINT CHR$(186) 

.OCATE 3, 1 
'RINT CHR$(200) ; 

•OR i% - 1 TO 78 
PRINT CHR$ ( 205 ) ; 
fEXT i% 

'RINT CHR$ ( 188 ) 

TestDataFileError is used to trap user input file name errors . 

)N ERROR GOTO TestDataFileError 
.OCATE 2, 3 
•ow% = 2 
:ol% = 24 

:olor 11, 0 

’RINT "Test Data File Name? 

:olor 12, 0 

MLITestDataFile$ contains the test data in a format compa table with that g 
INPUT MLITestDataFile$ 

)PEN MLITestDataFile$ FOR INPUT AS #1 
)N ERROR GOTO 0 

:olor 9, 0 
.OCATE 3, 1 
’RINT CHR$ ( 204 ) 

.OCATE 3, 40 
’RINT CHR$ ( 203 ) 

.OCATE 3, 80 
’RINT CHR$ ( 185 ) 

.OCATE 4, 1 
PRINT CHR$ ( 186 ) 

'.OCATE 4, 40 
PRINT CHR$( 186 ) 

.OCATE 4, 80 
PRINT CHR$ ( 186 ) 

.OCATE 5, 1 
PRINT CHR$ ( 200 ) 

.OCATE 5, 40 
PRINT CHR$ ( 202 ) 

.OCATE 5, 80 
PRINT CHR$ ( 188 ) 

POR i% = 2 TO 79 
LOCATE 5, i% 

IF i% <> 40 THEN PRINT CHR$(205) 

■TEXT i% 

r 

[nputBumperElasticModulus : 

.OCATE 4, 2 

:olor 11, 0 

PRINT "Bumper Elastic Modulus (MPa)? 

:olor 12, 0 

INPUT " " , BumperElasticModulus# 

IF BumperElasticModulus# <= 0# THEN 
LOCATE 6, 1 
COLOR 11, 9 

PRINT "Sorry - Bumper Elastic Modulus Must Be > Zero I" 

GOTO InputBumperElasticModulus 
2ND IF 

' 175 

InputBumperMassDensity: ORIGIN A I 

L0CATE 4 ' 42 OF 
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COLOR 11, 0 

PRINT " Bumper Mass Density (kg/nT3)? 

COLOR 12, 0 

INPUT BumperMass Density# 

IF BumperMa s s Dens i ty # <= 0# THEN 
LOCATE 6, 1 
COLOR 11, 9 

PRINT "Sorry - Bumper Mass Density Must Be > Zero l 
GOTO Input BumperMass Density 
END IF 

BumperSoundSpeed# = SQR(BumperElasticModulus# / BumperMassDensity# ) 

9 

' numdata% stores the number of data records in the database. 

numdata% = 0 

' The following variables store averages of the database records . 

BumpThkAve# = 0 

BumpStandOf f Ave# “ 0 

PressWallThkAve# “ 0 

ProjDiaAve# = 0 

ImpAngAve# = 0 

ProjVelAve# = 0 

BumpMaj Axis Ave# = 0 

BumpMinAxisAve# = 0 

MLIHoleDiamAve# = 0 

MLIMassLossAve# = 0 

PressWallMa jAxisAve# = 0 

PressWallMinAxisAve# = 0 

r 

'Modify top line of boarder. 

COLOR 9, 0 
LOCATE 5, 1 
PRINT CHR$ ( 204 ) 

LOCATE 5, 80 
PRINT CHR$ (185) 

9 

' Scroll through data. 

VIEW PRINT 6 TO 12 

9 

DO WHILE NOT E0F(1) 

numdata% = numdata% + 1 
INPUT #1, Dummy $ 

INPUT #1, TestID$ ( numdata% ) 

INPUT #1, DataSource$(numdata%) 

INPUT #1, Tes tDate$ ( numdata% ) 

COLOR 9, 0 
PRINT CHR$ ( 186 ) ; 

COLOR 11, 0 
PRINT ■ No . : " ; 

PRINT numdata%? 

PRINT " ID: 

PRINT TestID$ (numdata% ) ; 

PRINT " Source: " ; 

PRINT DataSource$(numdata%) ; 

PRINT " Date: 

PRINT TestDate$ (numdata% ) ; 

LOCATE CSRLIN, 80 
COLOR 9, 0 
PRINT CHR$ ( 186 ) 

INPUT #1, Dummy# 'Skip bumper material field. 

INPUT #1, BumperThickness#(numdata%) 
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BumpThkAve# - BumpThkAve# + BumperThickness#(numdata% ) 

INPUT #1, Bumpers tandOf f # ( numdata% ) 

BumpStandOf fAve# = BumpStandOf fAve# + BumperStandOff#(numdata%) 

INPUT #1, Dummy# 'Skip pressure wall material field. 

INPUT #1, PressureWallThickness#(numdata%) 

PressWallThkAve# = PressWallThkAve# + PressureWallThickness# (numdata% ) 
INPUT #1, Dummy# 'Skip projectile material field. 

INPUT #1, Pro jectileDiameter# (numdata% ) 

ProjDiaAve# = ProjDiaAve# + Pro jectileDiameter#(numdata% ) 

INPUT #1, ImpactAngle#(numdata%) 

ImpAngAve# = ImpAngAve# + ImpactAngle#(numdata%) 

INPUT #1, ProjectileVelocity#(numdata%) 

ProjVelAve# = ProjVelAve# + Pro jectileVelocity#(numdata% ) 

INPUT #1, BumperMa jorAxis# (numdata% ) 

BumpMa jAxisAve# = BumpMa jAxisAve# + BumperMa jorAxis# (numdata% ) 

INPUT #1, BumperMinorAxis# ( numdata% ) 

BumpMinAxisAve# = BumpMinAxisAve# + BumperMinorAxis#(numdata%) 

INPUT #1, MLIHoleDiam# ( numdata% ) 

MLIHoleDiamAve# = MLIHoleDiamAve# + MLIHoleDiam# ( numdata% ) 

INPUT #1, MLIMassLoss# ( numdata% ) 

MLIMassLossAve# = MLIMassLossAve# + MLIMassLoss# ( numdata% ) 

INPUT #1, PressWallMa jAxis#(nuradata% ) 

PressWallMa jAxisAve# = PressWallMa jAxisAve# + PressWallMa j Axis# (numdata% ) 
INPUT #1, PressWallMinAxis# (numdata% ) 

PressWallMinAxisAve# = PressWallMinAxisAve# + PressWallMinAxis# (numdata%) 
INPUT #1, Dummy $ 

LOOP 

f 

VIEW PRINT 

f 

' Build box for averages. 

LOCATE 12, 2 
FOR i% = 1 TO 78 
PRINT CHR$ ( 205 ) ; 

NEXT i% 

LOCATE 12, 1: PRINT CHR$(204) 

LOCATE 12, 40: PRINT CHR$(203) 

LOCATE 12, 80: PRINT CHR$(185) 

FOR i% = 13 TO 23 

LOCATE i% , 1: PRINT CHR$(186) 

LOCATE i% , 40: PRINT CHR$(186) 

LOCATE i% , 80: PRINT CHR$(186) 

NEXT i% 

LOCATE 24, 1: PRINT CHR$(200); 

FOR i% « 1 TO 78 

IF i% <> 39 THEN PRINT CHR$(205); 

IF i% = 39 THEN PRINT CHR$(202); 

NEXT i% 

PRINT CHR$( 188 ) 

' Calculate and print out parameter averages . 

BumpThkAve# = BumpThkAve # / numdata% 

LOCATE 13, 3 
COLOR 11, 0 

PRINT "Ave. Bumper Thk (in):"; 

COLOR 15, 0 

PRINT USING "##.###'' ~ ~ ~ " ; BumpThkAve#; 

BumpStandOf fAve# = BumpStandOf fAve# / numdata% 

LOCATE 13, 42 
COLOR 11, 0 

PRINT "Ave. Bump. Stand Off (in):"; w/ 
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COLOR 15, 0 

PRINT USING "##.### AAAA "; BumpS tandOf f Ave# ; 
PressWallThkAve# = PressWallThkAve# / numdata% 

LOCATE 15, 3 
COLOR 11, 0 

PRINT "Ave. Pres Wall Thk (in):"; 

COLOR 15, 0 

PRINT USING "## . ### AAAA " ; PressWallThkAve#; 

ProjDiaAve# = ProjDiaAve# / numdata% 

LOCATE 15, 42 
COLOR 11, 0 

PRINT "Ave. Proj. Dia. (in):"; 

COLOR 15, 0 

PRINT USING "##.### AAAA "; ProjDiaAve#; 

ImpAngAve# ■ ImpAngAve# / numdata% 

LOCATE 17, 3 
COLOR 11, 0 

PRINT "Ave. Impact Angle (deg):"; 

COLOR 15, 0 

PRINT USING "##.### AAAA "; ImpAngAve#; 

ProjVelAve# = ProjVelAve# / numdata% 

LOCATE 17, 42 
COLOR 11, 0 

PRINT "Ave. Proj. Vel. (km/sec):"; 

COLOR 15, 0 

PRINT USING "##.### "; ProjVelAve#; 

BumpMaj Axis Ave# = BumpMa jAxisAve# / nuradata% 

LOCATE 19, 3 
COLOR 11, 0 

PRINT "Ave. Maj. Bumper Hole (in):”; 

COLOR 15, 0 

PRINT USING "##.### AAAA "; BumpMa jAxisAve#; 
BumpMinAxisAve# = BumpMinAxisAve# / numdata% 

LOCATE 19, 42 
COLOR 11, 0 

PRINT "Ave . Min . Bumper Hole ( in ) : " ; 

COLOR 15, 0 

PRINT USING "##.### AAAA "; BumpMinAxisAve#; 
MLIHoleDiamAve# = MLIHoleDiamAve# / numdata% 

LOCATE 21, 3 
COLOR 11, 0 

PRINT "Ave. MLI Hole Diam. (in):"; 

COLOR 15, 0 

PRINT USING "##.### AAAA " ; MLIHoleDiamAve#; 
MLIMassLossAve# = MLIMassLossAve# / numdata% 

LOCATE 21, 42 
COLOR 11, 0 

PRINT "Ave. MLI Mass Loss (grams):"; 

COLOR 15, 0 

PRINT USING "##.### "; MLIMassLossAve#; 

PressWallMa jAxisAve# = PressWallMa jAxisAve# / numdata% 
LOCATE 23, 3 
COLOR 11, 0 

PRINT "Ave. Maj. P.Wall Hole (in):"; 

COLOR 15, 0 

PRINT USING "## .### AAAA " ; PressWallMa jAxisAve#; 
PressWallMinAxisAve# = PressWallMinAxisAve# / numdata% 
LOCATE 23, 42 
COLOR 11, 0 

PRINT "Ave. Min. P.Wall Hole (in):"; 
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DLOR 15, 0 

RINT USING "##.### PressWallMinAxisAve#; 

DCATE 25, 28 
DLOR 12, 0 

RINT "press any key to continue"; 

D 

DOP WHILE INKEY$ = "" 'Press any key to continue 

RuleOpt calculates the function coefficients and stores them in vector a#(l t 
Promt user for optimizer parameters . 

VLL OptParameters(InputIteration#, InputAlpha#) 

RuleOpt uses a modified Powell's method (as developed by W.K.Rule) to optimal 
\LL RuleOpt ( Inputlteration# , InputAlpha#, RSquaredValues#( ) ) 

ShowCoef ficients displays the calculated coefficients. 

■ILL ShowCoef ficients (RSquaredValues# ( ) ) 

Request the user for parameter values to be used for damage predictions . 
redictValue : 

DLOR 11, 0 
LS 

DCATE 1, 1 

RINT "ENTER DATA FOR DESIRED PREDICTION: " 

nputBumperThickness : 

DCATE 5, 1 
DLOR 11, 0 

RINT "Bumper Thickness (in): 

DLOR 12, 0 

NPUT "", PredictBumpThick# 

F PredictBumpThick# <= 0# THEN 
LOCATE 6, 1 
COLOR 11, 9 

PRINT "Sorry - Bumper Thickness Must Be > Zero I" 

GOTO InputBumperThickness 
ND IF 

DCATE 6, 1: PRINT " 

nputBumperStandOf f : 

DCATE 7, 1 
DLOR 11, 0 

RINT "Bumper Stand-Off (in): 

DLOR 12, 0 

NPUT "", PredictBumpStandOf f # 

F PredictBumpStandOf f# <= 0# THEN 
LOCATE 8, 1 
COLOR 11, 9 

PRINT "Sorry - Bumper Stand-Off Must Be > Zero!" 

GOTO InputBumperStandOf f 
ND IF 

DCATE 8, 1: PRINT " 

iputPressureWallThickness : 

DCATE 9, 1 
DLOR 11, 0 

RINT "Pressure Wall Thickness (in): 

DLOR 12, 0 

NPUT "", PredictPressWallThick# 

F PredictPressWallThick# <= 0# THEN 
LOCATE 10, 1 
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COLOR 11, 9 

PRINT "Sorry - Pressure Wall Thickness Must Be > Zero l" 

GOTO InputPressureWallThickness 
END IF 

LOCATE 10, 1: PRINT " 

I 

InputPro jectileDiameter: 

LOCATE 11, 1 
COLOR 11, 0 

PRINT "Projectile Diameter (in): 

COLOR 12, 0 

INPUT " ", PredictPro jDia# 

IF PredictPro jDia# <= 0# THEN 
LOCATE 12, 1 
COLOR 11, 9 

PRINT "Sorry - Projectile Diameter Must Be > Zerol" 

GOTO InputPro jectileDiameter 
END IF 

LOCATE 12, 1: PRINT " 

t 

InputlmpactAngle : 

LOCATE 13, 1 
COLOR 11, 0 

PRINT "Impact Angle (degrees): 

COLOR 12, 0 

INPUT " " , Predict ImpAngle# 

IF PredictlmpAngle# < 0# THEN 
LOCATE 14, 1 
COLOR 11, 9 

PRINT "Sorry - Impact Angle Must Be >= Zero!" 

GOTO InputlmpactAngle 
END IF 

LOCATE 14, 1: PRINT " 

/ 

InputPro jectileVelocity: 

LOCATE 15, 1 
COLOR 11, 0 

PRINT "Projectile Velocity (km/sec): "? 

COLOR 12, 0 

INPUT "", PredictPro jVel# 

IF PredictPro jVel# <= 0# THEN 
LOCATE 16, 1 
COLOR 11, 9 

PRINT "Sorry - Projectile Velocity Must Be > Zero!" 

GOTO InputPro jectileVelocity 
END IF 

LOCATE 16, 1: PRINT " 

t 

LOCATE 25, 28 
COLOR 12, 0 

PRINT "press any key to continue"; , 

DO 

LOOP WHILE INKEY$ ■ "" 'Press any key to continue 

f 

' MakePrediction evaluates the damage functions at user input values. 

CALL MakePrediction 

r 

' Show predictions on the screen, and associated function R-squared values ( 
CLS 
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OLOR 11, 0 

RINT "Calculated Results* 

OLOR 10, 0 

RINT "(Function R-Squared Values Given In Brackets)" 

If the impact angle for the prediction is zero, then set the minimum 
and maximum bumper hole diameters to be equal to their average . 

F Predict ImpAngle# = 0# THEN 

AverageDiameter# = ( BumpMaxDiaPred# + BumpMinDiaPred# ) / 2# 

BumpMinDiaPred# = AverageDiameter# 

BumpMaxDiaPred# = AverageDiameter# 

ND IF 

OCATE 5, 1 
OLOR 12, 0 

RINT "Minimum Bumper Hole Diameter (in): "; 

OLOR 15, 0 

RINT USING ”####.#### "; BumpMinDiaPred#; 

OLOR 10, 0 

RINT USING " (##.####) RSquaredValues# ( 1 ) 

OCATE 7, 1 
OLOR 12, 0 

RINT "Maximum Bumper Hole Diameter ( in ) s " ; 

OLOR 15, 0 

RINT USING "####.#### BumpMaxDiaPred#; 

OLOR 10, 0 

RINT USING " (##.####) "; RSquaredValues# ( 2 ) 

.OCATE 9, 1 
OLOR 12, 0 

'RINT "MLI Hole Diameter (in): " ; 

OLOR 15, 0 

'RINT USING "####.#### "; MLIDiaPred#; 

OLOR 10, 0 

•RINT USING " (##.####) "; RSquaredValues# ( 3 ) 

.OCATE 11, 1 
OLOR 12, 0 

’RINT "Pressure Wall Hole Diameter (in): " ; 

OLOR 15, 0 

•RINT USING "####.#### "; PressureWallDiaPred# ; 

OLOR 10, 0 

'RINT USING " (##.####) "; RSquaredValues# ( 4 ) 

Allow the user to make multiple predictions from the same set of coefficie 
.OCATE 24, 1 
OLOR 11, 0 

'RINT "Do You Wish To Enter Data For Another Prediction (y/n)? "; 

OLOR 12, 0 
NPUT " " , Answer$ 

OLOR 15, 0 

F Answer$ = "" THEN Answer$ = "y" 
jiswerS = LCASE$( Answer $) 

F Answer$ = "y" THEN GOTO PredictValue 
!ND 

This subroutine traps database input file errors. 

'estDataFileError : 

OLOR 12, 0 
.OCATE 4, 1 

•RINT "Please Re-enter File Name (or enter QUIT to stop)" 
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LOCATE row% , col% 

PRINT " 

LOCATE row%, col% 

INPUT MLITestDataFile$ 

MLITestDataFile$ ■ UCASE$ (MLITestDataFile$ ) 

IF MLITestDataFile$ - "QUIT" THEN END 
LOCATE 4, 1 
PRINT " 

RESUME 
QuitRunning : 

END 

RETURN 

SUB DisplayConvergence (col%, rsquared#) 

' This subroutine displays the effectiveness of the coefficient optimizer fo 

' This subroutine is designed to fill the screen with data, rather than scro 

' Only rows 3 thru 23 are used to display the data. 

CurRow% = CSRLIN 
CurCol% - POS(O) 

IF CurRow% - 24 THEN 
CurRow% = 5 
col% = col% + 10 
IF col% >80 THEN 
col% = 1 
COLOR 15, 0 

Here we clear the screen of data if it is full. 

FOR clrline% = 5 TO 23 

LOCATE clrline%. Is PRINT " 

NEXT clrline% 

COLOR 15, 9 
END IF 
END IF 

' rsquared# is the current value of the coefficient of determination (R~2) o 

LOCATE CurRow%, col% 

PRINT USING " #.#### rsquared# 

ENP SUB 

SUB MakePrediction 

' This subroutine uses the calculated function coefficients a#() to make pre 

t 

• Convert impact angle to radians . 

angle# = Predict ImpAngle# * pi# / 180# 

r 

' BumpMinDiaPred# is the predicted value of the bumper minimum hole diameter 

BumpMinDiaPred# = a#(l) * (PredictProjVel# / Bumper SoundSpeed# ) * a#(2) 
BumpMinDiaPred# = BumpMinDiaPred# * (PredictBumpThick# / PredictProjDia#) * a# (3 
BumpMinDiaPred# = BumpMinDiaPred# * ( COS ( angle# ) ) A a# ( 4 ) + a# ( 5 ) 

BumpMinDiaPred# * BumpMinDiaPred# * PredictProjDia# 

IF BumpMinDiaPred# < 0# THEN BumpMinDiaPred# - 0# 

t 

' BumpMaxDiaPred# is the predicted value of the bumper maximum hole diameter 

BumpMaxDiaPred# = a# ( 6 ) * (PredictProjVel# / BumperSoundSpeed# ) A a#(7) 
BumpMaxDiaPred# = BumpMaxDiaPred# * (PredictBumpThick# / PredictProjDia#) A a# (8 
BumpMaxDiaPred# = BumpMaxDiaPred# * ( COS ( angle# ) ) A a#(9) + a# (10) 
BumpMaxDiaPred# = BumpMaxDiaPred# * PredictProjDia# 

IF BumpMaxDiaPred# < 0# THEN BumpMaxDiaPred# = 0# 

t 

' MLIDiaPred# is the predicted value of the MLI hole diameter. 

MLIDiaPred# = a# (11) * (PredictProjVel# / BumperSoundSpeed#) ~ a# (12) 

MLIDiaPred# = MLIDiaPred# * (PredictBumpThick# / PredictProjDia#) * a# (13) 
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IDiaPred# = MLIDiaPred# * (PredictBumpStandOf f# / PredictPro jDia# ) A a#(14) 
IDiaPred# = MLIDiaPred# * ( COS ( angle# ) ) " a# (15) + a# (16) 

IDiaPred# = MLIDiaPred# * PredictPro jDia# 

MLIDiaPred# < 0# THEN MLIDiaPred# = 0# 

PressureWallDiaPred# is the predicted value of the average pressure wall h 
essureWallDiaPred# = a# (17) * (PredictPro jVel# / BumperSoundSpeed# ) * a# (18) 
essureWal IDiaPred# = PressureWallDiaPred# * (PredictBumpThick# / PredictPro jDi 
essureWallDiaPred# = PressureWallDiaPred# * (PredictBumpStandOf f# / PredictPro 
essureWallDiaPred# = PressureWallDiaPred# * (PredictPressWallThick# / PredictP 
essureWallDiaPred# = PressureWallDiaPred# * ( COS ( angle# ) ) ~ a#(22) + a#(23) 
essureWallDiaPred# = PressureWallDiaPred# * PredictPro jDia# 
PressureWallDiaPred# < 0# THEN PressureWallDiaPred# = 0# 

D SUB 

B ObjectiveFunction (0310%, code%, xtry#(), objective#, rsquared#) 

This subroutine calculates the objective function for the optimizer. 

Here the objective function is the coefficient of determination (R'‘2) of t 
objective# is the objective function. 

xtry#() is a vector of trial coefficient values used by the optimizer. 

:alc% = 0 means do not calculate R , '2, calc% = 1 means do calculate R~2. 
ijective# = 0# 

avemeasured# is the average measured value of the dependent parameter (use 
-emeasured# = 0# 

numdata% is the total number of records in the database. 

>R datacount% = 1 TO numdata% 

Convert impact angle to radians . 
angle# = ImpactAngle# (datacount% ) * pi# / 180# 

code% equal to 1 means treat bumper hole minor diameter function. 

IF code% = 1 THEN 

measured# is the measured value of the dependent variable, 
measured# = BumperMinorAxis# (datacount% ) / ProjectileDiameter# (datacount% ) 
calculated# is the calculated value of the dependent variable, 
calculated# = xtry#(l) * (ProjectileVelocity#(datacount%) / BumperSoundSpe 
calculated# = calculated# * (BumperThickness#(datacount% ) / Pro jectileDiam 
calculated# = calculated# * ( COS ( angle# ) ) ~ xtry#(4) + xtry#(5) 

END IF 

code% equal to 2 means treat bumper hole major diameter function. 

IF code% = 2 THEN 

measured# = BumperMa jorAxis# (datacount% ) / ProjectileDiameter# (datacount% ) 
calculated# = xtry#(l) * (ProjectileVelocity#(datacount%) / BumperSoundSpe 
calculated# = calculated# * (BumperThickness#(datacount%) / Pro jectileDiam 
calculated# = calculated# * ( COS ( angle# ) ) ~ xtry#(4) + xtry#(5) 

END IF 

code% equal to 3 means treat MLI hole diameter function. 

IF code% = 3 THEN 

measured# = MLIHoleDiam# (datacount% ) / Pro jectileDiameter# (datacount% ) 
calculated# = xtry#(l) * (ProjectileVelocity#(datacount%) / BumperSoundSpe 
calculated# = calculated# * (BumperThickness#(datacount%) / Pro jectileDiam 
calculated# = calculated# * ( Bumpers tandOf f#(datacount% ) / Pro jectileDiame 
calculated# = calculated# * (COS (angle# ) ) A xtry#(5) + xtry#(6) 

END IF 

code% equal to 4 means treat pressure wall average hole diameter function. 
IF code% = 4 THEN 

AverageDiameter# = (PressWallMinAxis#(datacount%) + PressWallMa jAxis#(data 
measured# = AverageDiameter# / Pro jectileDiameter# ( da tacount% ) 
calculated# = xtry#(l) * (Pro jectileVelocity#(datacount% ) / BumperSoundSpe 
calculated# = calculated# * (BumperThickness#(datacount% ) / Pro jectileDiam 
calculated# = calculated# * (BumperStandOff #(datacount% ) / Pro jectileDiame 
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calculated# - calculated# * (PressureWallThickness# (datacount% ) / Project 
calculated# ■ calculated# * ( COS ( angle# ) ) A xtry#(6) + xtry#(7) 

END IF 

objective# = objective# + (measured# - calculated#) ~ 2 
avemeasured# * avemeasure# + measured# 

NEXT da ta count % 

t 

IF calc% = 1 THEN 

avemeasured# = avemeasured# / numdata% 
value# = 0# 

' Here value# is determined which is used in the calculation of R~2. 

FOR datacount% = 1 TO numdata% 

IF code% = 1 THEN 

measured# * BumperMinorAxis#(datacount%) / ProjectileDiameter#(datacou) 
END IF 

IF code% = 2 THEN 

measured# = BumperMajorAxis#(datacount%) / ProjectileDiameter# (datacoui 
END IF 

IF code% = 3 THEN 

measured# = MLIHoleDiam#(datacount%) / ProjectileDiameter# (datacount%) 
END IF 

IF code% ■ 4 THEN 

AverageDiameter# = (PressWallMinAxis#(datacount%) + PressWallMa jAxis#(c 
measured# = AverageDiameter# / ProjectileDiameter# ( da tacount%) 

END IF 

value# = value# + (measured# - avemeasured#) A 2 
NEXT datacount% 

rsquared# = 1# - objective# / value# 

END IF 
END SUB 

SUB OptParameters ( Inputlteration#, InputAlpha#) 

COLOR 11, 0 
CLS 

LOCATE 1, 1 
PRINT "Enter the 
COLOR 14, 1 

PRINT "ITERATION PARAMETER"; 

COLOR 11, 0 

PRINT " for the function coefficient optimizer. " 

PRINT "Values in the range " ; 

COLOR 14, 1 

PRINT " (10 to 1000) "? 

COLOR 11, 0 

PRINT " are acceptable, 

COLOR 14, 1 

PRINT " 20 is recommended." 

COLOR 11, 0 

PRINT "High values will tend to produce better results but longer execution" 
PRINT "times." 

PRINT 

PRINT 

COLOR 14, 1 

INPUT " ITERATION PARAMETER? ", Inputlteration# 

f 

COLOR 11, 0 
LOCATE 10, 1 
PRINT "Enter the 
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COLOR 11, 0 

PRINT "for the function coefficient optimizer." 

PRINT "Values in the range " ; 

COLOR 14, 1 
PRINT " (0.1 to 3) 

COLOR 11, 0 

PRINT " are acceptable, "; 

COLOR 14, 1 

PRINT " 1 is recommended." 

COLOR 11, 0 

PRINT "High values will tend to reduce the chance of getting trapped in a" 

PRINT "local minimum (rather than the global minimum) but will tend to reduce" 
PRINT "the chance of precisely locating the global minimum." 

, PRINT 
PRINT 

COLOR 14, 1 

INPUT ” SEARCH DOMAIN PARAMETER? ", InputAlpha# 

LOCATE 25, 28 
COLOR 12, 0 

PRINT "press any key to continue"; 

DO 

LOOP WHILE INKEY$ = "" 'Press any key to continue 

END SUB 

SUB RuleOpt (Inputlteration#, InputAlpha#, RSquaredValues# ( ) ) 

' This subroutine finds optimal values for the prediction function coefficii 

' Optimal in the sense that function R A 2 values are minimized (nonlinear le; 

' The optimization technique is based on a modified Powell's method and is 

' described in the following reports 

' Rule, W.K., "ROCOPT - A User Friendly Interactive Code to Optimize 

' Rocket Structural Components," NASA CR-183837, 1989. 

/ 

DIM search# (1 TO 7, 1 TO 7 ) 'Matrix of columns which are search vectors 
DIM searchnew#(l TO 7) 'New search vector generated as a vector sum of previous 
DIM alpha# (1 TO 7) 'Search vector multiplier 
DIM x#(l TO 7) 'Design variables 

DIM xtry#(l TO 7) 'Trial values of design variables to check if objective funct 
seed% - ((TIMER * 65536) / 86400) - 32768 'Seed% is the seed number of the rand 
RANDOMIZE seed% 

ON KEY ( 1 ) GOSUB QuitRunning 
KEY ( 1 ) ON 

' If code%=l then find coefficients for BumpMinDiaPred# function. 

' If code%=2 then find coefficients for BumpMaxDiaPred# function. 

' If code%=3 then find coefficients for MLIDiaPred# function. 

' If code%=4 then find coefficients for PressureWallDiaPred# function. 

code% = 1 


FindCoef f icients s 

searchsequence% = 1 'This keeps track of the number of seaches run for a given 
Findcoef f icients 1 : 
iteration# = Inputlteration# 
alphatry# = InputAlpha# 


COLOR 15, 0 
CLS 

COLOR 11, 1 

IF code* = 1 THEN PRINT "Coefficient of 

IF code% = 2 THEN PRINT "Coefficient of 

IF code% = 3 THEN PRINT "Coefficient of 

IF code% = 4 THEN PRINT "Coefficient of 

LOCATE 3, Is COLOR 0, 15 


Determination for Bumper Hole Minimum D 
Determination for Bumper Hole Maximum D 
Determination for MLI Hole Diameter Fun 
Determination for Pressure Wall Hole Di 
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PRINT "CONDUCTING FUNCTION COEFFICIENT SEARCH SEQUENCE* 
COLOR 14, 0 

PRINT USING " # "; searchsequence% ; 

COLOR 0, 15 
PRINT " " ; 

COLOR 12, 0 

PRINT " PRESS FI TO QUIT 
COLOR 0, 15 
PRINT " 

COLOR 15, 0 

PRINT , , 

col% = l'This variable keeps track of the column position. 

COLOR 15, 9 


IF code% = 1 THEN numvar% = 5 

IF code% = 2 THEN numvar% = 5 

IF code% = 3 THEN numvar% = 6 

IF code% = 4 THEN numvar% =7 

totals% = iteration# * numvar% 'totals% is the total number of search matrices t 
alphamult# - .01# A (1# / totals%) 'This factor is to reduce alphatry# 

' to 1/100 of it's initial value by the end of the iterations. 

' Initialize variables. 

IF searchsequence% = 1 THEN 
FOR i% = 1 TO numvar% 
x#(i%) = 0# 

NEXT i% 

END IF 

nums% * 0 'Counter for number of search matrices generated 
iteration% - 0 'Counter for number of search vectors used 
start s 

FOR it = 1 TO numvar% 'Generate the random search matrix 
FOR j % = 1 TO numvar% 

search# (i%, j%) = -1# + 2# * RND 
NEXT j% 

NEXT i% 

FOR it ■ 1 TO numvar% 'Normalize random search vectors to +/-1. 
smax# = ABS ( search# ( 1 , i%)) 

FOR j% - 2 TO numvar% , u 

IF ABS ( search# ( j%, i%)) > smax# THEN smax# - ABS ( search# ( j%, i%)) 

NEXT j% 

FOR jt = 1 TO numvar% 

search# (j%, i%) - search# (j%, i%) / smax# 

NEXT j% 

NEXT i% 


nexts s 

calc% = 0 'calc% ■ 0 means do not calculate R A 2, calc% = 1 means do calculate R A 
FOR it = 1 TO numvar% 

iteration% * iteration% +1 . 

FOR j% = 1 TO numvar% 'Check objective function a negative distance along th 
xtry# ( j % ) = x#( j%) - alphatry# * search#(j%, i%) 

NEXT j % 

CALL ObjectiveFunction(calc%, code%, xtry#(), objective#, rsquared#) 
Obackward# = objective# 

FOR jt ■ 1 TO numvar% 
xtry# ( j % ) « x#(j%) 

NEXT j % 

CALL Ob jectiveFunction ( calc% , code%, xtry#(), objective#, rsquared#) 

Ocurrent# ■ objective# 

FOR jt = 1 TO numvar% 'Check objective function a positive distance along sea 
xtry# ( j % ) - x#(j%) + alphatry# * search# (j%, i%) 
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NEXT j% 

CALL ObjectiveFunction(calc%, code%, xtry#(), objective#, rsquared#) 
Oforward# = objective# 

IF Ocurrent# < Obackward# AND Ocurrent# < Oforward# THEN 

alpha# (i%) = 0# 'Make no change if current position is better. 

GOTO next s 2 
END IF 

IF Obackward# >= Oforward# THEN 
alpha# (i%) => alphatry# 

ELSE 

alpha# (i%) = alphatry# * (-1) 

END IF 

FOR j% = 1 TO numvar% 

x#(j%) = x#(j%) + alpha# ( i% ) * search#(j%, i%) 

NEXT j% 

novt* q • 

searchnew#(i%) = 0# 

NEXT i% 

IF nuras% >= total s% THEN GOTO finish 
nums% = nums% + 1 

alphatry# = alphatry# * alphamult# 

FOR ■ 1 TO numvar% 
xtry# ( j% ) = x# ( j % ) 

NEXT j% 

C£llO% = 1 

CALL ObjectiveFunction(calc%, code%, xtry#(), objective#, rsquared#) 

CALL DisplayConvergence(col%, rsquared#) 
smax# = 0# 

FOR i% « 1 TO numvar% 'Generate the new search vector 
FOR j% = 1 TO numvar% 

searchnew# ( i% ) = searchnew#(i%) + search# (i%, j%) * alpha# (j%) 

NEXT j % 

IF ABS ( searchnew# ( i% ) ) > smax# THEN smax# = ABS ( searchnew# ( i %) ) 

NEXT i% 

testl# - (nums% * 1#) / (numvar% * 1#) 
test2% = INT((nums% * 1#) / (numvar% * 1#)) 
test3# = testl# - test2% 

IF test3# = 0 THEN 
GOTO start 
ELSE 

IF smax# = 0# THEN 'Regenerate search matrix if current one does no good 
GOTO start 
END IF 

FOR j% = 1 TO numvar% 

searchnew# ( j%) - searchnew# ( j%) / smax# 

NEXT j% 

FOR i% = 1 TO numvar% 

FOR j% = 1 TO numvar% 

IF i% < numvar% THEN search# (j%, i%) = search# (j%, i% + 1) 

IF i% = numvar% THEN search# (j%, i%) - searchnew# ( j % ) 

NEXT j% 

NEXT i% 

GOTO nexts 
END IF 
finish: 

IF code% - 1 THEN j% = 1 

IF code% - 2 THEN j% = 6 

IF code% - 3 THEN j% - 11 

IF code% - 4 THEN j% = 17 

FOR it » 1 TO numvar% 
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a#(j% + i% - 1) - x#(i%) 

NEXT i% 

searchsequence% - searchsequence% + 1 

IF sear chsequence % <= 3 THEN GOTO Findcoef f icientsl 


C€llc% * 1 

CALL Object iveFunct ion (calc %, code%, x#(), objective#, rsquared#) 

' RsquaredValues# ( ) stores the IT 2 values for each prediction function. 

RsquaredValues# (code%) ■ rsquared# 
code% — code% + 1 

IF code% <= 4 THEN GOTO FindCoef ficients 


KEY ( 1 ) OFF 
END SUB 


SUB ShowCoef ficients (RSquaredValues# ( ) ) 

COLOR 11, 0 

CLS 

LOCATE 1, 24 

PRINT "Calculated Function Coefficients" 


COLOR 12, 0 
LOCATE 3/1 

PRINT "Minimum Bumper Hole Diameter Coefficients; (R A 2 * 

PRINT USING "##.#### RsquaredValues# ( 1 ) ; 

PRINT ")" 

LOCATE 5, 1 

PRINT USING " ###•##### a#(l); a#(2); a#(3); a#(4); a#(5); 

t 

COLOR 13, 0 

LOCATE 9, 1 ... . „ /T> ^ 0 _ „ . 

PRINT "Maximum Bumper Hole Diameter Coefficients; (R z = , 

PRINT USING "##.#### RSquaredValues# ( 2 ) ; 

PRINT ")" 

LOCATE 11, 1 

PRINT USING " ###.##### a#(6); a#(7); a#(8); a#(9); a# (10); 

r 

COLOR 14, 0 

LOCATE 15, 1 „ 

PRINT "MLI Hole Diameter Coefficients; (R 2 - ; 

PRINT USING "##.#### RSquaredValues# ( 3 ) ; 

PRINT ")" 

LOCATE 17, 1 

PRINT USING " ###.##### a#(ll); a#(12); a#(13); a#(14); a#(15); a#(16) 

r 

COLOR 10, 0 
LOCATE 21 1 

PRINT M Pressure Wall Hole Average Diameter Coef ficients t (R 2 s / 

PRINT USING "##.#### " ; RsquaredValues# (4 ) ; 

PRINT ")" 

LOCATE 23, 1 

PRINT USING " #*#.##### a# (17); a# (18); a# (19); a# (20); a# (21); a#(22); 


LOCATE 25, 28 
COLOR 12, 0 

PRINT "press any key to continue" ; 

DO 

LOOP WHILE INKEY$ - "" 'Press any key to continue 


a#( 23 
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